THERMAL  EFEECTS    IN  STELLAR  PULSATIONS 


By 

WILLIAM  DEAN  PESNELL 


A  DISSERTATION  PRESENTED  TO  THE   GRADUATE   SCHOOL 

OF  THE   UNIVERSITY   OF   FLORIDA 

IN   PARTIAL    FULFILLMEOT   OF  THE    REQUIREMENTS 

FOR  THE    DEGREE    OF   DOCTOR   OF  PHILOSOPHY 


UNIVERSITY   OF   FLORIDA 
1983 


ACKNOWLEDGEMENTS 

I  would  like  to  thank  Dr.  J.  Robert  Buchler  for  his  assistance 
during  the  past  four  years  as  well  as  his  great  patience.  In  addition,  I 
would  like  to  thank  A.  N.  Cox  for  insight  into  the  physical  nature  of 
stellar  pulsations  and  Dr.  Oded  Regev  for  useful  suggestions  and 
discussions. 

Several  other  people  have  assisted  in  parts  of  the  work.  They 
include  Dr.  Robert  Coldwell,  who  assisted  in  the  writing  of  the  programs 
that  are  included  and  Sharon  Bullivant,  who  demonstrated  the  use  of  the 
word  processor.  I  would  also  like  to  thank  Marie  Jose  Goupil  for  reading 
the  rough  drafts  of  the  manuscript. 

Some  of  the  computations  in  this  work  were  performed  at  the  Los 
Alamos  National  Laboratory,  Los  Alamos,  New  Mexico,  while  I  was  a  Summer 
Graduate  Research  Assistant  there  in  the  summer  of  1982.  The  bulk  of  the 
computations  were  done  at  the  University  of  Florida,  primarily  under  the 
MUSIC  subsystem,  for  which  I  am  grateful. 


ii 


TABLE   OF   CONTENTS 

ACKNOWLEDGMENTS ii 

ABSTRACT iv 

CHAPTER 

I  INTRODUCTION 1 

I I  BACKGROUND 3 

III  THEORY 7 

IV  THE  THERMAL  OPERATOR 10 

Asymptotic   Expansions 16 

Asymptotic   Thermal   Eigenf unctions 22 

V  NUMERIC  RESULTS 33 

BW  Vulpeculae 39 

Cepheid  model 41 

VI  DISCUSSION 74 

VII  CONCLUS  IONS 91 

APPENDICES 

A  THE    ANALYTIC  EQUATION   OF  STATE 93 

B  THE    INITIAL  MODEL    INTEGRATOR 103 

BIBLIOGRAPHY 125 

BIOGRAPHICAL  SKETCH 129 


ill 


Abstract  of  Dissertation  Presented  to  the  Graduate  School 
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By 
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A  Sturm-Liouville  operator,  embedded  in  the  linearized  radial 

pulsation  equations,  is  used  to  analyze  the  thermal  effects  of  stellar 

pulsations.  Asymptotic  expansions  are  utilized  to  demonstrate  the  nature 

of  these  eigenf unctions  and  to  define  the  thermal  response  time  of  a 

stellar  model.  The  eigenvalues  and  eigenvectors  for  models  of  a  Cepheid 

and  Beta  Cephei  variable  are  given  showing  the  differences  in  the  two 

classes  of  variables.  The  time  scales  introduced  agree  with  earlier  ones 

in  their  common  use  of  predicting  the  location  of  possible  driving 

zones.  In  addition,  a  way  of  displaying  the  temperature  variations  in 

the  hydrogen  ionization  zone  without  the  customary  "spike"  is 

demonstrated,  and  a  new  class  of  pulsations,  called  "sudden"  modes,  is 

introduced. 
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CHAPTER  I 
INTRODUCTION 

The  use  of  time  scales  to  Interpret  phenomena  found  in  stars  has 
been  a  favorite  tool  in  astrophysics,  (see  Cox  and  Giuli  1968;  Saio  and 
Wheeler  1982;  and  Shibahashi  and  Osaki  1981).  Primarily,  they  are  used 
to  simplify  the  governing  equations  by  ignoring  time  dependences  that 
are  much  longer  or  averaging  those  that  are  much  shorter  than  the 
subject  under  study.  Examples  of  this  include  the  suppression  of 
hydrodynamic  and  possibly  thermal  phenomena  when  examining  nuclear 
evolution  ( Iben  1965),  nuclear  evolution  effects  in  stellar  pulsations 
(Christy  1966),  and  hydrodynamic  phenomena  when  studying  thermal  flashes 
(Henyey  and  Ulrich  1972).  These  time  scales  arise  by  recasting  the 
various  relationships  into  non-dimensional  forms  and  defining  the 
factors  multiplying  the  time  derivatives  as  position  dependent  time 
scales.  This  produces  a  good  approximation  for  the  discussion  of  some 
attributes  of  pulsations. 

However,  we  have  been  using  asjnnptotic  formulations  to  study 
stellar  pulsations  (Buchler  1978),  and  when  we  attempted  to  include  the 
outer  layers  (Buchler  and  Regev  1982a) ,  it  was  found  necessary  to 
understand  the  thermal  structure  in  a  more  quantitative  way.  Other 
authors  have  found  modes  that  are  dominated  by  thermal  effects  (Wood 
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1976;  King  1980;  and  Saio  and  Wheeler  1982)  which  have  been  called 
"strange  modes."  In  this  work,  a  mathematical  formalism  for  studying 
modes  with  these  characteristics  is  given. 

We  shall  clarify  what  the  thermal  response  of  a  star  is  and  show 
effects  in  pulsations  that  result.  A  Sturm-Liouville  operator,  embedded 
in  the  linear,  non-adiabatic  radial  pulsation  equations,  is  used  to  give 
what  will  be  called  the  thermal  eigenf unctions  of  a  star.  It  will  be 
shown  that  under  physically  realistic  assumptions,  the  eigenvalues  are 
real  and  represent  decays  back  to  the  initial  state,  here  assumed  to  be 
one  of  thermal  equilibrium. 

Asymptotic  forms  for  the  eigenvalues  and  eigenf unctions  are  shown, 
and  several  realistic  models  have  their  spectra  calculated.  A  discussion 
of  the  physical  requirements  for  the  quasi-adiabatic  approximation  for 
the  stability  coefficient  is  given,  along  with  a  way  of  evaluating  it 
without  the  introduction  of  cutoffs  in  the  integration  variable.  The 
lowest  order  pulsation  operator  is  shown  to  have  a  more  complicated 
structure  than  the  simple  adiabatic  operator,  when  the  outer  layers  are 
included.  This  operator,  and  the  accompanying  stability  coefficient,  are 
derived  for  a  general  star  and  compared  to  earlier  work.  A  new  class  of 
modes,  dominated  by  thermal  effects,  is  shown  to  exist. 


CHAPTER  II 
BACKGROUND 

The  theoretical  study  of  stellar  pulsations  can  be  divided  into 
several  general  areas,  first  is  the  adiabatic  wave  equation  and  the 
discovery  that  the  periods  predicted  by  this  theory  agree  with 
observations  (Ledoux  and  Walraven  1958;  Cox,  J, P.  1980).  Here  the  entire 
star  is  assumed  to  react  faster  via  an  acoustic  process  than  thermally, 
implying  that  stellar  pulsations  are  basically  sound  waves  propagating 
in  the  envelope.  In  much  of  a  star's  mass,  adiabaticity  is  a  very  good 
assumption  and  the  periods  calculated  will  be  correct  as  long  as  this  is 
true.  The  stability  analysis  of  an  adiabatic  mode  uses  the  quasi- 
adiabatic  approximation  (Ledoux  and  Walraven  1958;  Ledoux  1965), 
although  the  integrals  must  have  a  cutoff  introduced  to  prevent  the 
outer  layers  from  dominating  the  result.  These  calculations  demonstrate 
that  the  quasi-adiabatic  interior  would  damp  any  oscillations  unless  an 
active  source  of  driving  was  present.  Several  authors  (Cox  1960  and 
1958;  Zhevakin  1963  and  his  earlier  papers)  showed  that  the 
thermodynamic  work,  derivable  from  the  ionization  of  once  ionized  helium 
could  provide  the  necessary  driving  for  the  classical  Cepheid  variables. 
This  driving  mechanism  requires  that  some  of  the  star  be  non-adiabatic 
and  a  finely  tuned  relationship  between  the  location  in  the  star  of  the 
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ionization  zone  and  the  transition  from  adiabatic  to  non-adiabatic 
behavior  be  satisfied.  Predictions  of  this  theory  are  narrow  instability 
strips  in  the  H-R  diagram  (Cox,  J. P.  1980),  with  sharply  defined  blue 
edges  where  the  evolving  star  first  becomes  unstable. 

Other  classes  of  stars  are  driven  by  similar  mechanisms.  Both 
radial  and  non-radial  oscillations  in  the  ZZ  Ceti  variables  are  driven 
by  hydrogen  ionization  (Starrfield  et  al.  1982;  Winget,  Saio  and 
Robinson  1982).  Hot  white  dwarfs,  in  particular  PG1159-035,  have  been 
shown  to  be  unstable  to  radial  pulsations  in  Starrfield,  Cox  and  Hodson 
(1980),  and  this  analysis  has  been  extended  in  Starrfield  et  al.  (1983), 
to  Include  non-radial  motion  and  to  show  the  driving  is  due  to  the 
ionization  of  carbon  and  oxygen.  Variables  such  as  the  RR  Lyrae,  BL 
Herculae,  6  Scuti  and  W  Virginis  classes,  are  driven  primarily  by  helium 
ionization  although  hydrogen  does  participate  in  the  destabilization 
(Cox,  J. P.  1980). 

Once  the  destabilizing  mechanism  was  known,  computer  codes  were 
developed  that  followed  the  nonlinear  behavior  of  pulsating  stars 
(Christy  1964;  King  et  al.  1964),  and  reproduced  many  observed  features 
of  the  Cepheid  and  RR  Lyrae  variables.  Details  such  as  the  asymmetry  of 
the  light  and  velocity  curves  and  the  phase  shift  of  the  luminosity 
maximum  with  respect  to  the  radius  minimum  were  explained  using  large 
(and  expensive)  codes.  A  number  of  new  problems  arose  from  this  work, 
two  major  ones  being  the  mass  discrepancy  in  Cepheids  (Cox,  A.  N.  1980a) 
and  the  lack  of  driving  in  Beta  Osphei  variables  (Osaki  1982;  Lesh  and 
Aizenman  1978).  In  the  attempt  to  understand  these  difficulties,  fully 
non-adiabatic,  but  linear,  stability  analyses  were  developed  (Baker  and 
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Kippenhahn  1962;  Castor  1971;  and  Saio,  Winget,  and  Robinson  1983). 
These  analyses  have  the  advantages  of  being  much  easier  to  use  and  less 
expensive  than  the  nonlinear  codes.  They  are  easily  modified  to  include 
more  physics,  such  as  convection  (see,  for  example.  Baker  and  Gough 
1979). 

The  linear  analysis  shows  only  whether  a  given  model  is  self- 
exciting  or  damping.  An  unstable  model  in  the  linear  regime  is  a 
candidate  for  further  study  in  a  hydrodynamic  code.  However,  the  use  of 
the  linear,  non-adiabatic  eigenvectors  in  a  perturbation  expansion  is 
limited  due  to  their  non-periodic  nature.  To  understand  how  non- 
linearities  affect  an  oscillation,  without  using  a  hydrodynamic  computer 
code,  a  purely  real  spectrum  of  eigenvalues  is  required  in  first  order. 
The  adiabatic  eigenf requencies  satisfy  this  requirement,  but  treat  the 
outer  layers  in  an  extremely  poor  fashion. 

As  the  bulk  of  the  star  is  quasi-adiabatic,  Buchler  (1978) 
developed  a  formalism  using  two-time  asymptotic  techniques  modeled  after 
Cole  (1968)  and  Kuzmak  (1959).  This  method  links  the  linear  analysis 
with  the  nonlinear  calculations  by  providing  a  consistent  expansion  in 
powers  of  the  amplitude  (Buchler  and  Perdang  1979).  The  requirement  of 
periodicity  is  inherent  in  an  expansion  of  this  type,  but  the  payoff  is 
the  capability  to  study  a  pulsation  amplitude  as  it  starts  at  an  initial 
value  and  evolves  to  a  final  state.  Additional  problems  that  can  be 
studied  include  modal  selection,  mode  switching,  and  the  feedback  of  the 
pulsation  on  the  thermal  evolution  of  the  star,  all  of  which  are  open 
questions  in  this  field  (Regev  and  Buchler  1981;  Simon,  Cox  and  Hodson 
1980). 
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In  the  simplest  version  of  the  theory,  the  quasi-adiabatic  results 
are  recovered  in  the  zero  amplitude  limit,  but  at  the  cost  of 
introducing  a  cutoff  in  the  integrals.  When  this  was  done,  it  became 
obvious  that  these  layers  were  being  incorrectly  treated  and  that  a  more 
reliable  method  needed  to  be  developed.  A  straightforward  way  of 
including  these  layers  was  presented  in  Buchler  and  Regev  (1982a),  where 
the  structure  of  the  eigenvectors  was  determined  by  the  condition  of 
constant  luminosity  above  the  cutoff  point  in  the  model.  While  giving 
good  results  for  stars  of  the  Beta  Cephei  class  (Pesnell,  Regev  and 
Buchler  1982),  cooler  stars,  like  Cepheids ,  were  not  treated  much  better 
than  with  the  simpler  adiabatic  approximation.  The  present  work  was 
undertaken  to  quantify  the  thermal  behavior  of  pulsations.  The  results 
presented  are  encouraging  and  explain  or  clarify  several  phenomena  found 
in  stellar  pulsations. 


CHAPTER  III 
THEORY 

The  structure  and  evolution  of  stars  is  governed  by  a  set  of  four 
equations:  the  conservation  of  momentum,  energy,  and  mass,  and  a  flux 
equation.  Other  relationships  are  constituitive ,  such  as  the  equation  of 
state,  opacity  (see  Appendix  A),  and  nuclear  energy  generation  rates. 
Written  in  a  Lagrangian  form  which  uses  the  interior  mass  as  the 
independent  variable,  these  equations  are  (see  Clayton  1968,  ch.  6; 
Appendix  B  below) 


^=   4up  (3-1) 

0  m 


dr  Gm  /^^P-zon  ft   o^ 

= 4Tir<^  -r —  =  g(r,S)  (3-2) 

o  o  0  m 

d    t^  r^ 


T^.<i(p,T)-^-=IHr.S)  (3-3) 


«»)  -  -  (*-^)^  f  H'  "-*> 


The  symbols  are  defined  here  to  be 

m  =  mass  beneath  the  point  of  interest 
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r  =  radius  from  the  center  of  the  star  to  this  point 

p  =  density 

T  =  temperature 

P  =  pressure,  P  =  P(p,T) 

S  =  entropy,  S  =  S(p,T) 

L  =  luminous  flux 

K  =  opacity,  k  =  k(p,T) 

q  =  nuclear  energy  generation  rate,  q  =  q(p,T). 

We  have  assumed  radiative  heat  transfer  in  the  diffusion 
approximation  and  that  all  constitutive  quantities  are  given  by  analytic 
expressions. 

To  study  this  system,  the  equations  are  solved  with  all  time 
derivatives  set  to  zero  to  find  an  initial  model  in  hydrostatic  and 
thermal  equilibrium.  The  equations  are  then  linearized  in  r  and  S  and 
combined  to  give  two  operator  equations  with  two  unknowns,  6r  and  6S. 
Written  in  operator  form  these  are 


77=  (^•^^-(||-^^^-A-5r-B.6S        (3-5) 


d  6S  _/5  k\.^^  +(14). 6S  =  C.6r  +  D.6S.       (3-6) 


d  t    \d   r/     V5  S 


It  is  well  known  that,  with  the  appropriate  boundary  conditions, 
the  A  operator  is  of  the  Sturm-Liouville  class  and  yields  a  symmetric 
matrix  operator  (Castor  1971;  Buchler  and  Regev  1982b). 
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The  operator  A  is  the  adiabatic  pulsation  operator  (6S=0)  which  can 
be  written  as  (see  Cox,  J.  P,  1980) 


dl  l^,fp''[^^  II  *  I'^d?  [(3^i-*>^]  -  ^^'  ^  1  ^ 


dt^' 


and  in  the  Schrodinger  form  of  this  equation  (Morse  and  Feshbach  1953, 
ch.  6) ,  the  new  independent  variable  has  the  dimensions  of  time  and 
represents  the  acoustic  travel  time  from  a  point  in  the  star  to  the 
surface  and  is  given  by 


t     (r)  =   /  di 
acous      r*^ 


r^(r)  P(r) 


p(r) 


=  ^/  dr/c(r)  (3-7) 


where  c(r)  is  the  local  adiabatic  sound  speed  and  R^   is  the  equilibrium 

radius  of  the  star.  The  fundamental  radial  pulsation  period  is  roughly 

given  by  t     (R    ),  the  sound  travel  time  from  the  core  to  the 
acous   core 

surface.  As  the  acoustic  time  corresponds  to  the  time  for  a  dynamic 

response  to  occur,  t     (r)  can  be  identfied  as  the  dynamic  time  scale 
'   acous  ■^ 

as  a  function  of  position  in  a  star. 

The  nice  properties  of  A  come  from  the  fact  that  it  is  of  second 
order  in  the  spatial  derivatives  (the  time  dependence  is  assumed  to  be 
exp(  itot  )),  Looking  at  the  D  operator,  the  same  second  order 
derivatives  are  present.  We  now  show  that  there  is  a  representation 
where  D  has  the  same  properties  as  A,  and  a  natural  definition  of  the 
thermal  time  scale  will  result. 


CHAPTER  IV 
THE  THERMAL  OPERATOR 

For  an  optically  thick  medium,  the  thermal  operator  D  can  be 
written  in  a  continuum  form  as 


^  d   6S  ,  „      .  6S  d 


d   L 


AC  u    1    tA-I   d   6S/C 

T  /-    \  /■  /  \    ^S    ,    ^  ,    .  /d    lnT\  V 

L(m)(4   -   K„)  -pr-  +   Uta) 


eq       T   C^^    6q  y  d  m  y    dm 

(4-1) 


Here  q  =(^  "^|  ,  -3 — -^   =  q(  P   ,T   ),  i.e.  an  initial  model  in  thermal 
^T  \dlnTy  '  d  m     ^  ^eq  eq 

balance,  k     =  L  "!1  j  ,   C  ={-r-^)  »  and  all  derivatives  are  evaluated  at 

the  equilibrium  densities  and  temperature.  This  form  for  D  has  been 

linearized  at  constant  density  (radius)  and  is  strictly  valid  at  points 

deep  in  the  star  so  that  t;  >  1. 

We  assume  a  time  dependence  of  the  form  exp(  -at  ) ,  and  note  that 

the  operator  is  of  second  order  in  the  spatial  derivative.  Therefore,  an 

integrating  factor,  \x   ,  exists  that  transforms  D  into  a  self-adjoint 

operator.  This  factor  can  be  written  as 

r°,,      s  d  InT  , 
„  J  (4  -  K   )  dm 

m^"^        T   d  m 

li(m)  =  e   "  (4-2) 


where  ra^  is  a  reference  mass,  here  chosen  to  be  m    ,  the  mass  of  the 
0  core 
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core  where  nuclear  burning  occurs.  Defining  a  new  eigenvector  as 

d  In  L 

I   =  6S/C  and  \  =  — ; ^^ ,  equation  4-1  becomes 

^  V  dm 

d    K_  ,  An    A-l 


p  c„.  c  -  [x(4  -.,-,,)- 7^  k  -  3^  [  .(VS'J- 


*---  -  dm. 

eq      ' 


(4-3) 


Appropriate   boundary   conditions    (which  make   4-3   self-adjoint)    are 


and 


i-^^  =  a   ^(m,)  +  b  4^(m^)   =   0  (4-4a) 

dm  *  dm     " 


C(mQ)   =   0.  (^-^b) 


The  surface  boundary  condition  (4-4a)  is  equivalent  to  the  assumption 
that  at  the  surface  there  is  no  radiation  incident  from  the  outside 
(Castor  1971).  In  this  form,  4-4a  shows  that  the  thermal  readjustment 
time  of  the  outer  surface  is  infinitely  short.  As  the  weight  function  on 
the  left-hand  side  is  positive  definite,  and  the  function  multiplying 
the  derivative  on  the  right-hand  side  is  negative  definite,  4-4  is  of 
the  Sturm- Li ouville  class  of  operators  (Boyce  and  DiPrima  1969,  ch.  11). 

Comparing  4-3  to  the  standard  Sturm- Liouville  equation,  since 

A-1 

>  0  and  uTC  /L   >  0  ,  all  eigenvalues  are  real  and 
V   en 


/d  InfT 


v   eq 
positive.  With  each  eigenvalue  a     is  associated  an  eigenvector  ^^  whose 

number  of  nodes  corresponds  to  the  numeric  position  of  the  {a^}  in  an 

increasing  ordering  of  the  eigenvalues  starting  with  a^.    In  addition. 
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the  eigenvectors  are  non-degenerate,  orthogonal  with  respect  to  the 
weight  function  jiC  T/  L  ^  ,  and  can  be  normalized  so  that  (  m^  is  the 
total  mass  of  the  star) 


0 


/    fi  T  C  I.    l.    dm  =   6.  ., 


eq 


Due  to  their  Sturm-Liouviile  character,  these  eigenvectors  form  a 
complete  set  for  describing  the  6S  readjustment  in  a  star.  In 
combination  with  the  adiabatic  radial  wave  functions  (Cox,  J.  P.  1980), 
a  complete  set  for  describing  spherically  symmetric  disturbances  of  a 
star  is  found,  provided  that  equations  3-1  through  3-4  are  always  valid. 

To  understand  the  structure  of  the  eigenvectors,  the  approximation 

of  neglecting  nuclear  burning  can  be  made.  This  reduces  the  equation  to 

a  version  suitable  for  discussing  thermal  effects  in  a  stellar  envelope. 

Without  nuclear  burning  the  initial  luminosity  is  constant  and  the 

analysis  is  halted  at  a  mass  (  m     )  outside  the  region  where  this 
■'  core  ^ 

assumption  is  not  valid.  As  the  structure  of  the  ^  's  will  show,  this  is 
equivalent  to  assuming  we  are  looking  at  "short"  thermal  time  scales 
that  do  not  probe  the  core.  Radial  stellar  pulsations  are  an  envelope 
phenomenon,  so  this  assumption  is  justified. 

Neglecting  nuclear  burning  and  defining  L^  to  be  the  static 
luminosity  entering  the  envelope  from  below,  equation  4-3  becomes 


L^        vn       dm   n   dm 


I  d  m  J  dm 


(4-5) 


Following  Ledoux  and  Pekeris  (1941),  we  multiply  eq.  4-5  by  C  and 
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integrate  over  the  interval  m^=  m    <  m  <  m . .  A  self-adioint  system 
*'  0   core       *  J      ^ 

satisfies  a  minimum  principle  and  any  well  behaved  eigenvector  can  be 
used  to  estimate  the  lowest  eigenvalue,  a„.  Choosing  Cr)(™)  -   1   for  all 


m,  this  estimate  is 


m^    d  <^  m^  d  K^ 

J  i^^"")  TIT-   dm         „  J  ^i  T^r-  dm 


^        dm             m^''    d  m 
o„  =  -  Lj.  : =  -  L^  


r  *  <U> 

/  H  T  C  dm  ^^      ,,  ,- 

m  •"  *^    v  (4-6) 


If  the  thermal  time  is  defined  as  the  reciprocal  of  a„  ,  this  gives  the 
overall  thermal  time  of  the  envelope  as 


_        <U> f,    -,. 

t^,  = r .  (4-7) 

th      ^     m^   d  K^ 

*         I      a  —: dm 

m  -^      dm 

There  is  an  apparent  mathematical  inconsistency  with  this 

expression.  When  the  opacity  derivative,  <„,  is  constant  in  the 

envelope,  the  thermal  time  is  undefined.  However,  comparing  4-6  to  the 

expression  for  the  adiabatic  radial  pulsation  fundamental  frequency 

(Ledoux  1965), 


0  ^  ''^^^  dT  t  (  3  r^-  4  )  P  ]  dm 


o/ 


r'^  dm 


d  P 
they  have  the  same  form  only  if  — - —  =  0  ,  so  that  the  formula 

d  P         "^  ^1 
for  t^,  is  incomplete.  If  - —  =  0  and  - —  =  0,  then  cjo_=  0,  but  the 
th         '^  dm         dm  0 

variations  of  P  are  large  in  a  star,  varying  many  orders  of  magnitude, 
while  r  can  vary  only  from  1  to  5/3  and  the  neglect  of  the  pressure 
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derivative  would  not  be  correct.  The  inclusion  of  nuclear  burning  in  4-7 
remedies  this  situation  and  the  net  thermal  time  for  the  star  is  then 


th.nuc         m^  d  k 

L*  0^   I't^^^  -  K^-   q^)  -  ^— -  ]  dm 

which  allows  for  the  luminosity  variations  and  shows  that  the  apparent 

infinite  value  of  t^,  from  4-7  when  k_  is  constant  in  the  star  is  not 
th  T 

real.  The  terra  allowing  for  nuclear  burning  will  always  give  a  finite 

value  for  t  .      .   When  k_  is  constant  in  the  envelope,  the  assumed 
th,nuc  T 

form  of  C  is  a  poor  approximation  and  gives  an  unrealistic  result.  The 

positivity  of  t  ,      is  difficult  to  guarantee,  but  a  simple  case  would 
'^  ■'  th.nuc  J 

d  K 

require  that  4  -  k„  -  q„  >  0,  — ; >  0,  and  —r —  <  0  ,  all  of  which  can 

T    T       dm  dm 

be  satisfied  in  most  of  a  star's  mass,  although  in  some  nuclear  burning 
regions,  q  can  be  much  greater  than  zero,  and  the  opacity  variations 
are  discussed  below.  The  inclusion  of  nuclear  burning  as  it  relates  to 
stability  of  stellar  models  is  discussed  in  Hansen  (1978)  and  Gabriel 
(1972).  In  this  work,  this  area  will  not  be  discussed  as  the 
simplifications  that  result  from  the  neglect  of  the  core  make  the 
interpretation  of  the  results  easier. 

The  opacity  is  of  the  Kramers'  form  (  k  ~-3.5  <  0  )  near  the 
surface  but  inside  the  ionization  zones.  Near  the  core,  electron 
scattering  will  dominate,  (  k  — >0),  so  that  k„  decreases  going 
outward  in  the  star.  Therefore,  except  in  a  small  region  near  the 
ionization  zones,  t  ,   will  have  only  positive  contributions  in  the 
integral.  The  mass  of  the  ionization  zones  is  small  compared  to  the 
total  mass  of  a  star,  so  that  the  negative  contributions  to  the  integral 
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will  not  cause  a  negative  overall  time  scale.  In  addition,  each  point  in 
the  star  is  weighted  by  the  ratio 


/    s    ~   I    T(,m;  \     T 
|a(m)   '   '  '  » 


T( 


which  gives  the  interior  regions  far  more  importance  than  the  outer 
layers,  further  reducing  the  influence  of  the  negative  portions  of  the 
integrand.  Physically  t   represents  the  time  necessary  to  remove  the 
|J.  weighted  internal  energy  <U>  by  the  radiative  luminosity  present, 
allowing  for  the  position  (and  therefore  temperature)  dependence  of  the 
opacity  derivative. 

Another  thermal  time  can  be  found  from  equation  3-3  as  the  time 
necessary  to  remove  the  material  thermal  energy  at  the  static  luminosity 
present  in  the  envelope.  Writing  this  in  a  integral  form  (Saio,  Winget 
and  Robinson  1983)  yields 


^h  n=  J  *T  Cdm/L..  (4-8) 


'th,p   m 


Comparing  the  two  definitions,  there  are  several  differences.  The  value 
t^j^   is  counting  only  the  thermal  energy  in  the  matter  that  makes  up 

the  star.  In  equation  4-7,  the  formula  emphasizes  the  radiation 

4 
(  IJ.  ~  T  )  and  the  method  of  energy  transport,  as  well  as  giving  the 

interior  stellar  regions  more  weight  than  the  outer  layers.  These 

differences  are  minor  and  are  overshadowed  by  a  major  shortcoming 

of  t^^.  The  definition  in  4-8  can  be  generalized  to  calculate  a  thermal 

time  as  a  function  of  mass  by  replacing  the  lower  limit  in  the  integral 
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by  the  desired  mass.  Values  for  t  ,    found  are  always  positive  and  can 
■^  th,p 

be  interpreted  as  time  scales.  However,  when  the  same  technique  is  done 

with  t  ,  ,  the  value  returned  is  not  positive  definite  due  to  the 
th 

ionization  zones  near  the  surface.  Mathematically,  this  generalization 
means  that  CCni')  =  1>  m'>m  and  C(ni' )  -   0»  ni'<m,  a  Heaviside  step 
function  at  m.  The  derivatives  of  C  in  4-5  will  yield  delta  functions  at 
m,  adding  contributions  to  t  .  that  are  not  easily  understood.  To 
alleviate  these  difficulties,  t  .  shall  be  redefined  in  terms  of  the 
Independent  variable  that  results  from  analyzing  eq.  4-3  by  standard 
asymptotic  methods. 

Asymptotic  Expansions 

Following  the  technique  outlined  in  Itorse  and  Feshbach  (1953,  ch. 
6)  we  make  the  substitutions: 


m    TC   ,  T  ^  1/ 
1    r  r     V  d  InT  f/o  . 

^  =  T   J  L  ~  T"  ~~A J   dm  , 

J  m;;  "-   L.   dm-" 


0 

m.     TC 

/   [ ^ 

"»0      li*    dm 


,=    j"*[_i^  V^t/2,.  , 


1  2    ,2 
k.  =  J  a  , 


and 
w(x) 


k^!Mki''.iiJi  th 


t-<¥Sf-v] 


1  "   L  2j^d  m/   V  ^    J^   -^   -  "T 

d  x^ 


rl    .  1/.       A   .2  ,2  TC   d  m 
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to  transform  eq.  4-3  to 

,2 

I    +  (  k2  -  w(x)  )  y  =  0.  (4-9) 

d  X 

1/- 
The  units  of  J  are  (time)^  ,  so  a  natural  definition  of  the  overall 

thermal  time  of  a  stellar  envelope  is 

t^h=  J2.  (4-10) 

We  note  that  the  integrating  factor  [i   does  not  appear  in  this 
definition.  As  ^  is  an  artificial  quantity  introduced  for  mathematical 
convenience,  this  is  physically  appealing. 

From  the  equilibrium  model  (B-4  in  Appendix  B  below)  there  is  a 
relationship  between  L^  and  dlnT/dm.  Inserting  this  into  eq.  4-10  gives 


m^   3  C  K       i- 


=  [   /   ( )'2dm  ]2  (4-lla) 


^^  '   ""O    (47ir2)2  ac  T 


or,  using  the  continuity  equation  (3-1),  R^  =  r(m_)  and  R^^  =  r(m^). 


K^   J  '-   ^  1/ 

^h  =  t  J  ( — ^-—y^p^-  ]'-  (^-iib) 

0    a  c  T^ 


This  definition  shows  the  effect  of  both  the  internal  energy  that  must 
be  dissipated  and  the  method  of  dissipation.  Rewriting  4-ll(b)  in  a 
slightly  different  way 


R.    3  K  T  C    1/ 
0     c  aT^ 
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or,    defining   U  =   T   C   , 

'  °     gas  V 


R^     U 


"0       "rad 


^h  ■  i^y<'^^>'''^^'         <*-'» 


we  have  a  interesting  definition  for  t  .  .  The  photon  diffusion  length  is 

present  (  l/<p  ),  so  increasing  the  optical  thickness  increases  the 

relaxation  time.  The  3/c  term  can  be  interpreted  as  the  diffusion 

velocity  of  the  photon;  however,  this  is  dicussed  below.  From  the  ratio 

of  internal  energies,  it  is  seen  that  in  the  central  regions  where  this 

ratio  is  large,  the  contributions  to  t   increase.  Near  the  surface, 

where  the  density  and  opacity  approach  zero  and  the  temperature  becomes 

constant,  the  radiation  dominates  and  the  integrand  becomes  small. 

If  the  approximation  of  a  monatomic  ideal  gas  with  black  body 

radiation  is  made,  a  simple  relationship  exists  for  the  ratio 

U   /  U   .  Defining  n  to  be  the  reciprocal  of  the  mean  molecular  weight 
gas   rad 

(or  the  molar  density  in  moles/gr) ,  the  pressure  of  the  gas  can  be 
written  as 


a  4 
P  =  P   +  P   ^  =  nRpT  +  ^  T, 
gas   rad  J 


or,  introducing  P  ( Chandrasekhar  1967, ch.  6), 


P    =  BP  and  P   ,  =  (l-p)P. 
gas  rad 
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The  internal  energy  of  the  gas  can  be  written  as 


3         4 
U  =  U   ,  +  U    =  T  nRT  +  aT  /p. 


'as    rad   2 


so  that 


U    =  -^  P   /p  =  4  P  P/P 
gas   2  gas  "^   1   '^ 


and 


"rad  =  3^ad/P  =  (1-P)P/P- 


With  these  assumptions,  the  ratio  becomes 


U   /U   ,  =  1/2  P/  ( 1  -  p)  . 
gas   rad    ^  ^  ^ 


This  mixture  has  a  thermal  time  given  by 


R  1 

^h  =  If  U^  *  t  <P  P/(1-P)r2dr  }2.  (4-13) 


If  P  is  constant  in  the  envelope 


'th  -|j-T4l[af*<'^''^-']^-  <*-'^> 
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Here  we  see  that  the  thermal  time  is  a  diffusion  length 


^*    1/ 

1  =  [  u/   (<P)  2dr  ]2  (4-15) 

^0 


divided  by  a  diffusion  velocity 


V  =  c  I  (1  -  p)/p,  (4-16) 


where  c  is  the  speed  of  light  in  vacuo. 

The  diffusion  velocity  must  be  less  than  c;  therefore,  the  value 

of  p  must  be  greater  than  0,4.  This  formula  shows  that  in  stars  where 

P  is  of  order  1,  the  thermal  times  are  large  as  it  takes  a  long  time  for 

a  photon  to  move  through  the  material.  In  a  star  that  is  radiation 

dominated,  the  thermal  time  is  short  as  the  photons  can  move  freely 

through  the  envelope.  These  interpretations  are  qualitative,  depending 

on  the  assumed  nature  of  U   .  They  are  suitable  for  estimating  whether 

gas 

a  star  undergoes  adiabatic  pulsations  or  moves  in  some  way  that  shows 
large  thermal  effects.  A  star  with  a  very  tenuous  envelope,  and 
therefore  small  p,  would  be  expected  to  show  large  thermal  effects. 
Stars  of  this  type  include  the  R  Corona  Borealis  class  (King  1980;  Saio 
and  Wheeler  1982),  as  well  as  a   Cygnus  (Lucy  1976;  Pesnell  1981). 

Returning  to  the  original  form  of  t   ,  a  position  dependent  time 
scale  can  be  defined  by  replacing  the  lower  limit  by  the  desired  radius. 
The  resultant  formula  represents  the  thermal  travel  time  from  that  point 
to  the  surface  and  can  be  written  as 
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R^   3  K  C 


t,u(r)  =  [  /  \ ^/^2p  dr  ]2.  (4-17) 

r     a  c  T^ 


In  figure  4-1,  we  compare  this  thermal  time  with  t  ,   ,  also 
^  '^  th,p' 

generalized  to  a  position  dependent  form  as  discussed  above  (Cox,  J. P. 

1980;  Saio,  Winget,  and  Robinson  1983).  The  time  scales  are  plotted 

versus  the  logarithm  of  temperature  (with  the  surface  of  the  star  at  the 

left  of  the  diagram)  for  a  model  of  BW  Vulpeculae  and  in  figure  4-2,  for 

a  generic  Oepheid,  discussed  in  chapter  5.  Except  for  the  constant 

separation,  it  is  readily  seen  that  the  two  functions  agree  quite  well 

in  the  interior  but  diverge  from  one  another  at  the  surface.  This  effect 

is  a  consequence  of  t     containing  only  the  material  thermal  energy, 
tn ,  p 

which  for  a  hot  star  is  quite  considerable  and  the  drop  to  zero  at  the 

surface  is  due  only  to  the  integral  formulation  for  t^,   .  In  the 

th,p 

equation  for  t   ,  the  method  of  energy  transport  reduces  to  the  grey 

atmosphere  approximation  at  the  surface,  and  the  thermal  time  must 

therefore  approach  zero  more  smoothly  for  any  class  of  star. 

This  disagreement  does  not  change  the  normal  use  of  the  thermal 

time.  Several  classes  of  variable  stars  are  driven  by  ionization  zones, 

including  Gepheids  (Cox,  J. P.  1980)  and  ZZ  Ceti  (Starrfield  et  al.  1982; 

Winget,  Saio,  and  Robinson  1982).  For  this  driving  mechanism  to 

function,  the  ionization  zone  must  coincide  with  the  point  in  the  star 

where  the  transition  from  (quasi-)adiabatic  (constant  entropy)  motion  to 

highly  non-adiabatic  (constant  luminosity)  motion  occurs  (Cox  and 

Whitney  1958).  This  transition  zone  is  found  at  the  radius,  R   , 

where  t  ,   (R   )  ~  11^,  the  thermal  time  is  of  the  same  order  as  the 
th,p   tr     0' 
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period  of  the  pulsation.  The  expression  for  t   is  derived  from  a 

frequency,  and  the  location  of  this  region  would  be  where 

a{R  )  =  t   (R   )   ~  cjj  =  2Ti/n  ,  which  is  analagous  to  critical  damping. 

The  constant  separation  in  the  figures  is  approximately  In;  so  both 
formulations  predict  the  same  location  for  the  transition  zone,  provided 
this  region  is  below  the  atmosphere. 


Asymptotic  Thermal  Eigenfunctions 

Once  the  normalized  coordinate  has  been  defined,  the  asymptotic 
forms  for  the  eigenvalues  and  eigenfunctions  of  eq.  4-9  can  be 
calculated.  Here,  asymptotic  implies  that  k^  =  J^  a  >w(x),  for  all 
values  of  x  considered.  Due  to  the  surface  effects,  this  is  not  a  good 
assumption,  but  an  important  feature  in  the  eigenfunctions  can  be  seen 
in  this  approximation;  so,  it  is  discussed. 

The  asymptotic  form  of  eq.  4-9  is  the  harmonic  oscillator  and  has 
boundary  conditions 

a'y(l)  +  b'  ^(1)  =  0  (4-18a) 

d  x 

and 

y(0)  =  0  (4-18b) 

where  a'  and  b'  are  a  and  b  from  eq.  4-4  transformed  to  the  x  coordinate 
system.  Assuming  y(x)  =  A  sin(  kx  +  6  )  replaces  the  outer  boundary 
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condition  by  a  transcendental  equation  of  the  form 

tan(  k  +  6  )  =  b'/a'  k  (4-19) 

In  the  limit  k  >  1  ,  this  equation  has  roots  that  correspond  to  half- 
integer  multiples  of  n   (Butkov  1968,  ch.  9)  ,  k  =  (  -j  +  i  )  't. 
a  =  (  (  y  +  1  )    n   )2/j2,  and  we  shall  let  6  =  0  . 

The  physical  eigenvector  can  be  found  from  the  asymptotic 
expression  as 

5^(,,  .,_^Cx)/[-.^^(V?yf'*.  (4-20, 

In  Figure  4-3,  we  plot  the  x  variations  for  the  Beta  Cephei  model. 
For  log(T)<5.7,  the  shape  of  the  eigenvector  is  determined  solely  by  the 
denominator  in  4-20.  In  this  star  the  integrating  factor  decreases 
outward  and  the  net  effect  is  an  eigenvector  that  peaks  at  the  surface. 
(Figures  4-4  and  4-5).  The  minimum  in  \x   at  log(T)"'4.6  is  due  to  a 
maximum  in  k   at  that  temperature  caused  by  the  second  ionization  stage 
of  helium.  The  shoulder  in  y  at  this  temperature  is  also  from  this 
minimum. 

Moving  to  a  Cfepheid  model,  the  situation  changes  considerably.  In 
Figure  4-6,  the  thermal  coordinate  is  plotted,  and  we  see  that  x~l 
until  log(T)  >4.5,  well  interior  to  the  hydrogen  ionization  zone.  When 
the  asymptotic  (n=10)  eigenvector  is  plotted  (Figure  4-7),  a  large  peak 
occurs  at  log(T)=3.9.  This  variation  cannot  be  due  to  the  cosine  term  in 
4-20  as  the  argument  is  essentially  1.0;  so,  this  feature  must  be  due  to 
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the  denominator,  whose  major  variations  come  from  the  integrating 

factor.  A  large  minimum  in  \i   can  be  seen  in  Figure  4-8,  where  the 

logarithm  of  |j.  is  plotted.  This  minimum  corresponds  to  the  maximum 

in  5  and  is  caused  by  the  large  positive  value  for  k     in  the  cool  side 

of  the  hydrogen  ionization  zone.  The  integrating  factor  has  a 

4-K^ 
temperature  dependence  of  roughly  T       ,  and  in  regions  where 

<  >  4,  \i   increases  outward  instead  of  decreasing  as  it  does  elsewhere. 

This  structure  corresponds  to  the  peak  in  the  temperature  eigenvector  in 

a  standard  linear,  non-adiabatic  pulsational  stability  analysis.  As  the 

numeric  results  will  demonstrate,  any  ionization  zone  will  cause  a  local 

peaking  of  the  eigenvector,  but  a  region  where  k  >  4  causes  an  effect 

that  is  best  removed  when  calculating  thermal  effects  in  these  stars, 

and  plotting  the  resultant  eigenvectors.  There  are  secondary  minima 

in  (1  at  log(T)~4.2  and  4.7  due  to  the  ionization  stages  of  helium.  Any 

effects  these  have  on  the  eigenvector  y  (Figure  4-7)  are  completely 

hidden  by  the  peak  from  the  hydrogen. 

However,  it  is  not  true  that  k^  =  J^o  >  w(x)   everywhere.  At  the 

surface,  T — >Tq  and  k  — >0,  so  that  w(x)  increases  outward.  While  the 

interpretation  of  the  "bumps"  in  the  eigenfunctions  as  maxima  in  ic 

will  be  valid  in  the  actual  numeric  results,  the  frequencies  and 

eigenfunctions  are  not  accurately  represented  in  these  expressions.  With 

this  knowledge  in  hand  we  now  move  on  to  solving  the  equations 

numerically  and  show  the  thermal  modes  for  several  stars. 
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CHAPTER  V 
NUMERIC  RESULTS 

In  this  section  the  eigenvalues  and  eigenvectors  for  two  stellar 
models  will  be  discussed.  In  each  case,  the  initial  model  was  found  by 
using  the  model  builder  described  in  Appendix  B,  with  mass,  luminosity, 
effective  temperature,  and  composition  appropriate  to  the  class  of  star 
the  model  represents.  One  model  is  of  BW  Vulpeculae,  a  hot  supergiant  of 
the  Beta  Cephei  class  of  variables,  and  the  other  is  a  generic  Cepheid, 

To  find  the  matrix  operator,  equations  3-3  and  3-4  are 
differentiated  at  constant  radius  (density)  and  written  in  a  form 
resembling  that  of  Castor  (1971): 


T  -^T^  =  AK2«(T6S).  (5-1) 

d  t 


Although  the  matrix  AK2  is  not  symmetric,  it  is  tridiagonal  and  a 
numeric  analog  of  the  integrating  factor  does  exist  (Wilkinson  1965), 
Applying  this  transformation  to  AK2  yields  the  D  operator  as  a  real, 
symmetric,  tridiagonal  matrix.  The  diagonal  transformation  matrix  is 
defined  as 


V9 
/2(1)  =  1 
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/2(i)  =  /2  (i_i)  [  AK2(I-1,3)/AK2(I,1)  P,    1=2, N, 

Vo 
and  the  matrix  D  is  found  by  using  |i  -^  as  a  similarity  transformation 

D  =  \i^  AK2  (i  2 

so  that  the  spectrum  of  D  is  identical  to  that  of  AK2.  In  a  convective 
region  of  a  star,  one  procedure  is  to  neglect  the  Lagrangian 
perturbation  of  the  convective  flux  (Saio,  Winget  and  Robinson  1983). 
This  will  have  no  direct  effect  on  the  definition  of    \i  "^    as  the  two 
matrix  elements,  AK.2(  1-1,3)  and  AK2(I,1),  are  multiplied  by  the  same 

fraction  of  the  luminosity  being  carried  by  the  radiation.  Any  gross 

V9 
features  of  [i '^  ,  such  as  peaks  in  ionization  regions,  will  still  be 

present  in  a  model  where  convection  is  permitted,  although  they  may  be 

extended  over  a  larger  region  due  to  the  smaller  temperature  gradient. 

The  transformation  has  one  side  benefit,  the  matrices  analogous  to  AG2 

and  AKl  in  Castor  (1971),  are  much  closer  in  magnitude  to  the  elements 

of  the  symmetric  matrices  A  and  D  and  have  smaller  variations  as  a 

function  of  position  in  the  model,  resulting  in  better  numeric  accuracy. 

Real  eigenvalues  are  guaranteed  by  the  symmetric,  tridiagonal  form 
of  the  matrix  operator  (Ralston  1967).  For  negative  definite  eigenvalues 
(representing  decays)  the  diagonal  elements  must  be  negative.  This  can 
be  seen  in  a  simple  example. 

A  negative  diagonal  element  means  that  increasing  the  entropy 
(temperature)  in  a  zone  increases  the  luminous  flux  out  of  the  zone 
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more  than  it  decreases  th»  flux  entering  it.  This  stabilizes  the  zone  as 
the  additional  energy  is  lost  and  the  temperature  decays  to  its  initial 
value.  A  simple  example  can  be  used  if  the  luminosity  gradient  had  a 
local  relationship.  Equation  5-1  would  then  be  written  as 


TT=-^^-^o> 


so  that 


S(t)  -  Sq  =  S^e^*^ 


If  a>0,  the  entropy  decays,  but  if  a<0  the  entropy  increases.  The  latter 
situation  is  possible  in  a  nuclear  burning  core  but  not  in  the  envelope. 
To  preserve  the  Sturm-Liouville  character  of  the  continuum  operator 
in  4-5  in  the  matrix  operator,  the  off-diagonal  elements  must  be 
positive.  This  condition  is  composed  of  two  separate  requirements.  The 
luminosity  exiting  a  zone  must  increase  if  the  temperature  of  the  zone 
above  is  decreased,  and  the  luminosity  entering  a  zone  must  also 
increase  if  the  temperature  of  the  zone  below  is  increased.  Both  of 
these  are  easily  satisfied  if  the  opacity  is  independent  of  temperature 
or  decreases  with  temperature,  (  k   <  0  ),  From  the  form  of  the 
continuum  operator  (4-5),  only  if  <  >  4  can  violations  of  these 
conditions  occur  and  these  regions  are  precisely  where  the  interpolation 
formulae  for  the  luminosity  are  essential.  In  the  initial  model  the 
temperature  ratio  is  not  allowed  to  exceed  a  given  number  in  the 
ionization  zones.  This  process  limits  the  variation  of  the  opacity  as 
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well  and  raises  the  upper  bound  for  positive  matrix  elements.  A  first 
order  expansion  of  the  super-diagonal  matrix  element  shows  that  for 


K^  <  l/ln(T(I-l)/T(I))~17  , 


these  matrix  elements  are  always  positive.  We  have  never  found  any 
negative  off-diagonal  matrix  elements  in  any  calculated  envelope  model. 
This  process  also  limits  the  variations  of  the  matrix  elements  so  that 
the  sum  of  the  off-diagonal  elements  is  less  than  the  diagonal  elements, 
permitting  the  use  of  Geschgorin's  Theorem  (Wilkinson  1965). 

The  boundary  conditions  are  incorporated  into  the  matrix  operator 
by  giving  certain  elements  of  the  matrix  special  functional  forms. 
Assuming  the  outer  surface  to  be  radiative  and  given  by  the  grey 
atmosphere  formula  is  equivalent  to  assuming  that  this  surface  has  an 
infinitely  short  thermal  response  time.  Therefore,  the  analytic  form  of 
this  boundary  condition  is  that  of  eq.  4-4,  even  though  the  diffusion 
approximation  is  not  very  accurate  here.  For  numeric  computations  this 
expression  is  unsuitable,  as  it  is  valid  at  a  point  exterior  to  the 
actual  surface  of  a  model,  and  we  have  assumed  that  the  exiting 
luminosity  is  given  by 


L(W-l)  =  4Tta  r2(N+1)  T^^^  ^  (.    x  +  -|  ) 


When  this  formula  is  linearized  and  inserted  into  the  operator  equation, 
the  outer  equation  reduces  to  one  of  first  order  and  closes  the  system. 
The  inner  boundary  condition  is  6S/C  — >  0  and  is  included  by  zeroing 
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any  elements  in  the  matrix  that  refer  to  5S/C  (0).  Again,  this  reduces 

the  appropriate  equation  (here  the  innermost  zone)  to  one  of  first 

order. 

As  the  thermal  time  of  the  outer  surface  is  zero,  and  the  overall 

thermal  time  (    ~  a^        )  is  normally  longer  than  a  pulsation  period,  the 

spectrum  of  D  must  contain  eigenvalues  which  bracket  the  radial 

fundamental  adiabatic  pulsation  frequency.  Here  is  the  information  that 

we  seek  regarding  the  transition  zone.  All  parts  of  a  star  that  can 

relax  faster  through  a  thermal  process  than  acoustically  will  move  in  a 

highly  non-adiabatic  fashion.  These  regions  of  a  star  will  be  called  the 

"sudden"  regions  and  correspond,  in  some  stars  (see  below),  to  the  outer 

layers  where  — r—^ =  0  in  a  fully  non-adiabatic  stability  analysis. 

dm 

The  inner  regions  of  a  star  have  thermal  response  times  longer  than  the 
fundamental  period  (  a  <  w„  )  and  are  referred  to  as  quasi-adiabatic. 

Eigenvalues  for  this  matrix  can  be  found  by  a  variety  of  means.  For 
the  low  order  modes  iterating  on  the  outer  boundary  condition  is 
possible  (Castor  1971),  but  the  structure  of  the  high  order  modes  is 
such  that  this  method  fails.  A  general  eigenvalue  routine  (Smith  et  al. 
1976)  is  used  to  find  all  eigenvalues  of  D  and  the  necessary 
eigenvectors  are  found  by  standard  Gaussian  back  substitution.  The 
eigenvectors  of  D  and  AK2  are  related  by  the  same  similarity 
transformation  as  the  matrices  and  the  more  desirable  set  can  be 
calculated,  depending  on  the  effective  temperature  (  T  ,-  ),  of  the  star 
under  study.  If  T  ^^  >  15,000  K,  then  the  hydrogen  ionization  zone  is 
absent  and  the  C  's  give  a  good  representation  of  the  thermal  modes. 
When  T  ,,  <  15,000  K,  this  ionization  zone  causes  a  spike  to  be  present 
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in  ^  as  shown  in  the  section  on  asyiiiptotic  eigenvectors  in  Chapter  4. 
To  remove  this  feature,  the  eigenvectors  of  D  are  calculated  and 
displayed.  This  results  in  a  much  clearer  picture  of  what  is  occurring 
in  the  surface  regions  of  the  star. 

We  cannot  look  at  the  long  decay  rate  modes  as  we  have  removed  the 
core  and  the  effect  it  has  on  the  eigenvectors.  The  lower  modes  can  have 
significant  amplitude  in  the  core,  but  the  higher  modes  are  effectively 
isolated  from  the  core  even  more  so  than  the  radial  pulsation  modes. 
This  causes  the  difficulty  noted  earlier,  the  number  of  nodes  in  a 
calculated  eigenvector  is  no  longer  an  indicator  of  it's  position  in  the 
spectrum;  the  size  of  the  numbers  causes  the  computer  to  underflow  and 
the  interior  portion  of  the  eigenvector  is  set  to  zero!  Physically,  this 
implies  that  the  exterior  can  react  to  perturbations  on  much  shorter 
time  scales  than  the  interior.  However,  the  outer  layers  are  linked 
quite  strongly  to  the  interior  through  the  low  order  modes.  When  a  heat 
pulse  is  inserted  at  the  base  of  the  envelope,  the  entire  model  is 
affected.  The  outer  layers  react  to  maintain  thermal  balance  with  the 
new  luminosity  from  the  heat  pulse.  The  eigenf unctions  of  the  diffusion 
equation  exhibit  an  infinite  propagation  velocity  for  perturbations 
(Morse  and  Feshbach  1953,  ch.  7),  causing  this  behavior.  The  asymptotic 
eigenfunctions  are  trigonometric  with  an  amplitude  that  increases 
outward  and  we  can  therefore  speak  of  the  thermal  time  of  a  point  in  the 
star  as  the  reciprocal  of  the  largest  eigenvalue  whose  eigenvector  has 
significant  amplitude  at  that  position.  This  may  seem  ambiguous,  but  as 
thermal  modes  can  be  used  to  expand  the  pulsation  equations  in  the 
thermal  times,  this  identification  will  enable  us  to  study  the  thermal 
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effects  in  an  oscillation.  In  particular,  the  quasi-adiabatic  stability 
coefficient  and  pulsations  that  are  dominated  by  thermal  effects  can  be 
defined. 

Concentrating  on  those  modes  with  have  eigenvalues  comparable  to 
that  of  the  radial  fundamental  pulsation  frequency,  we  now  examine  two 
stellar  models. 

BW  Vulpeculae 

The  physical  parameters  for  this  star  are  plotted  in  Figure  5-1  and 
the  eigenf unctions  6r/r,  6p/p  and  6L/L  for  the  radial,  adiabatic, 
fundamental  pulsation  mode  are  shown  in  Figure  5-2.  The  results  of  a 
linear,  non-adiabatic  radial  pulsation  analysis  (LNA)  are  shown  in  Table 
5-1.  To  find  a  realistic  model  for  this  star,  it  was  necessary  to 
integrate  into  the  nuclear  burning  core  slightly,  but  the  thermal  modes 

discussed  below  are  not  altered  by  this. 

-4    -1 
The  fundamental  frequency  ((jl)„)  is  3.9  *  10   sec   ,  and  the  four 

thermal  eigenvectors  with  eigenf requencies  closest  to  this  are  plotted 

in  Figure  5-3. 

Important  features  are  the  cutoff  near  log(T)=5.2  and  the  two 

peaks.  The  peak  near  log(T)=4.6  is  due  to  the  second  ionization  of 

helium  causing  a  local  maximum  in  k  .  At  log(T)~5.1,  a  secondary  peak 

associated  with  the  first  due  to  the  use  of  the  Rosseland  mean  opacity. 

The  weighting  function  in  this  mean  has  a  maximum  at  a  temperature 

roughly  four  times  that  where  the  ionization  occurs  and  causes  a  maximum 

in  K     at  this  temperature.  This  peak  has  been  studied  by  Stellingwerf 
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(1978)  as  well  as  Saio  and  Cox  (1980),  as  a  possible  driving  mechanism 

for  these  stars. 

The  luminosity  variations  (  6L/L  )  for  each  mode  are  shown  in 

Figure  5-4.  Near  the  surface  —r. =  0  ,  showing  the  domination  of  the 

dm 

short  thermal  times  in  this  region.  Moving  inward  (increasing  the 

temperature),  the  luminosity  variation  starts  to  decrease  at  log(T)=5.0, 

oscillates  due  to  the  nodal  structure  of  C  »  and  drops  to  zero  by 

log(T)=6.0.  These  modes  are  completely  isolated  from  the  nuclear  burning 

core  (log(T)>6.5)  as  stated  earlier. 

Figures  5-5  and  5-6  show  the  real  and  imaginary  parts  of  the 

entropy  variations  as  calculated  from  the  LNA  pulsational  stability 

analysis.  The  peaks  at  log(T)=4.6  and  5.1  as  well  as  a  sharp  cutoff 

above  log(T)=5.3  can  be  seen  in  both  components.  Each  feature  has  an 

analogy  in  the  thermal  eigenvectors  in  Figure  5-3.  The  modulus  of  the 

luminosity  variations  from  this  analysis  is  shown  in  Figure  5-7.  The 

change  from  — ; =  0  is  present,  albeit  at  a  slightly  higher 

dm 

temperature  than  in  Figure  5-4.  Some  of  this  is  due  to  contamination 
with  low  order  thermal  modes,  but  the  significant  amplitude  in  6L/L  (up 
to  log(T)=6.0)  is  due  to  the  adiabatic  density  fluctuations  in  this 
region. 

This  model  shows  that  the  shape  of  6S/C  in  an  LNA  pulsational 
stability  analysis  is  determined  by  the  thermal  modes  of  that  star.  As 
always,  the  model  is  stable  and,  at  present,  we  do  not  find  any  new  ways 
of  destabilizing  it.  The  structure  of  the  eigenvectors  serves  to 
demonstrate  and  clarify  the  physical  reasons  behind  a  transition  zone 
and  the  splitting  of  the  star  in  two  regions,  adiabatic  (log(T)>5.4)  and 
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instantaneous  thermal  balance  (or  sudden)  at  log(T)<5.4.  More 
interesting  behavior  can  be  seen  when  a  model  of  a  Cepheid  variable  is 
examined. 

Cepheid  model 

This  star  is  a  generic  Cepheid  whose  characteristics  were  obtained 
from  A.  N,  Cox  (1980a).  The  results  of  the  LNA  pulsational  analysis  and 
summary  of  the  input  data  are  found  in  Table  5-2.  A  graph  of  the 
physical  parameters  is  found  in  Figure  5-8  and  a  graph  of  the  adiabatic 
eigenvectors  6r/r,  6p/p  and  6L/L  is  shown  in  Figure  5-9.  The  model  is 
unstable  in  the  first  three  pulsation  modes  and  corresponds  to  a  model 
of  the  star  DL  Cassiopeiae.  The  results  presented  here  are  illustrated 
very  nicely  by  this  star  and  are  similar  to  those  of  other  Cfepheids 
calculated. 

In  Figure  5-10  and  5-11  the  four  thermal  modes  that  bracket  the 
pulsation  frequency  CjO_  ,  are  plotted.  These  are  the  eigenvectors  of  AK2 , 
and  we  see  almost  no  detail  (or  the  individual  modes)  due  to  the  spike 
in  the  hydrogen  ionization  zone.  The  origin  of  this  spike  was  discussed 
in  the  section  on  asymptotic  eigenvectors,  and  as  it  eliminates  any 
other  features  present,  we  choose  to  remove  it  from  the  eigenvectors. 
This  is  accomplished  by  using  the  eigenvectors  of  D,  as  shown  in  Figure 
5-12.  These  vectors  have  the  same  eigenf requencies  as  those  in  Figure  5- 
10,  but  much  more  detail  can  be  seen.  The  bumps  at  log(T)=3.82  and  4.12 
are  due  to  the  hydrogen  ionization  zone.  They  are  the  lower  slopes  of 
the  spike,  and  the  minimum  at  log(T)=4.0  is  a  result  of  modulation 
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effect  of  the  integrating  factor  being  a  maximum  at  this  point.  The 

helium  ionization  zone  is  visible  at  log(T)=4.6  as  an  inverted  peak,  but 

the  second  peak  at  log(T)=5.1  is  absent.  From  Figure  4-6,  where  the 

thermal  coordinate  for  this  star  is  plotted,  we  see  that  the  value  for  x 

starts  to  drop  away  from  one  at  this  point.  In  the  interpretation  of 

features  for  log(T)<4.7,  the  neglect  of  the  trigonometric  function 

in  ^  is  valid,  but  interior  to  this,  the  actual  nodal  structure 

of  t     will  start  to  dominate, 
n 

Examining  the  luminosity  variations  which  correspond  to  Figure  5-12 
(Figure  5-13),  the  driving  mechanism  responsible  for  destabilizing  this 
star  is  visible  (see  below).  To  accentuate  this,  the  variations 
for  a  <  0)   are  shown  in  Figure  5-14  and  5-15  and  a  >  co  in  Figure  5-16 
and  5-17.  In  the  low  frequency  results,  the  luminosity  looks  much  like 
that  of  BW  Vul  (Figure  5-4);  increasing  the  frequency  to  those  in  Figure 
5-16  and  5-17,  the  decrease  in  5L/L  near  log(T)=4.6  has  become  larger 
than  any  other  part  of  the  vector.  Here  is  the  physical  reason  for  the 
driving  zone  of  this  star. 

For  a  star  to  be  unstable  in  a  pulsation  mode,  a  region  where 

^  ^^^^  <  0  must  exist  (Cox,  J.  P.  1980).  This  is  seen  if  the  work 
d  m 

integral  is  writen  in  the  form  (Aizenman  and  Cox  1975;  Buchler  and  Regev 
1982b), 


=  nJ   (^3-  1)  (-^)   (-j-^)  dm  . 


The  quantity  T^-   1  is  always  positive  and  for  the  fundamental  mode 

— ^  <  0  (as  —  increases  outward)  (see  Figure  5-2);  thus,  when 
P  r 
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'^   >  0  the  star  is  stabilized  and  for  —3 <  0  the  star  is 

dm  dm 

unstable.  In  the  BW  Vul  model  (Figure  5-7)  — ^ —  is  zero  in  the 

atmosphere  and  6L/L  generally  decreases  with  increasing  temperature.  The 

thermal  modes  have  the  same  basic  luminosity  behavior  (Figure  5-4), 

showing  that  this  star  is  stable.  In  the  Cepheid  model,  the  pulsational 

6L/L  (Figure  5-26)  shows  a  definite  region  near  log(T)~4.6  where 

—   <  0  .  This  outward  decrease  of  6L/L  cannot  be  due  to  the 

d  m 

adiabatic  density  fluctuations  as  they  will  always  have  6L/L  increasing 

outward.  The  thermal  modes  show  a  change  in  structure  as  the  decay  rate 

is  incresased  through  the  value  a~(i).  The  change  in  the  shape  of  6L/L  in 

this  region  is  the  cause  of  the  destabilization  of  the  Cepheid. 

Calculations  of  the  BW  Vul  model  with  o/w  <  1  and  co/o  <  1  do  not  show 

the  same  change.  Examining  Figure  5-18  and  5-19  where  the  thermal  modes 

with  o/u)~10,  and  Figures  5-20  and  5-21  where  Wa~10,  we  do  not  see  the 

same  change  in  the  structure  of  the  thermal  modes.  This  leads  to  the 

conclusion  that  the  star's  thermal  structure  prevents  this  type  of 

driving  mechanism  from  functioning. 

The  location  of  the  transition  zone  in  each  model  should  be 

discussed  as  well.  The  BW  Vul  model  shows  a  definite  change  from 

adiabaticity  to  sudden  behavior  at  log(T)=5.2.  The  Cepheid  shows  a  less 

dramatic  change,  as  the  thermal  eigenvectors  do  not  show  an  extremely 

rapid  cutoff  at  the  transition.  However,  unlike  BW  Vul,  the  Cepheid  has 

a  region  where  —: <  0  that  begins  near  the  transition.  Also,  it 

dm 

would  seem  more  difficult  to  establish  a  transition  zone  for  this  model 

as  r- —  is  not  zero  until  well  above  the  helium  ionization.  In  fact, 

d  m 

the  luminosity  variations  from  the  pulsation  analysis  (Figure  5-26)  do 
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not  have  a  vanishing  derivative  until  above  the  hydrogen  ionization 
zone.  It  is  difficult  to  use  the  timescale  arguments  to  locate  the 
transition  zone  for  such  a  model.  The  use  of  a  «  a)_  to  partition  the 
thermal  modes  is,  in  this  case,  a  fairly  reliable  way  to  locate  this 
region.  For  a  <  a)_  ,  the  thermal  mode  has  a  6L/L  vector  that  is,  in 
general,  increasing  outward  and  has  only  a  small  dip  in  the  helium 
ionization  zone.  When  a  >  oo.,  6L/L  is  dominated  by  this  region, 
until  a     >  a)„  ,  when  the  amplitude  of  the  eigenvector  (5S/C  )  becomes 
small  enough  that  the  mode  is  isolated  from  that  region  of  the  star. 

To  see  the  benefits  of  the  integrating  factor,  the  real  and 
imaginary  parts  of  the  LNA  pulsational  analysis  are  shown  in  Figures  5- 
22,  5-23,  5-24  and  5-2  5.  In  Figures  5-22  and  5-23,  the  eigenvectors  of 
the  AK2  matrix  are  shown,  with  the  "spike"  due  to  the  hydrogen 
ionization  zone  dominating  the  graphs.  When  the  eigenvectors  of  D  are 
examined  in  Figures  5-24  and  5-25,  the  variations  have  been  greatly 
reduced.  In  5-22,  ( 6S/ C  )   ,  varies  by  two  orders  of  magnitude  and  the 
imaginary  part  changes  by  almost  three.  The  eigenvectors  of  D  vary  by 
less  than  one  order  of  magnitude  and  the  details  of  the  helium 
ionization  and  other  structures  can  be  seen. 

The  eigenvectors  shown  in  this  section  have  some  numeric  noise  that 
deserves  comment.  The  system  that  is  solved  here  is  identical  to  that 
used  in  the  LNA  stability  analysis,  and  any  noise  in  the  thermal  modes 
is  also  present  in  that  analysis.  The  "spikiness"  of  the  eigenvectors 
present  in  Figure  5-14  can  be  eliminated  by  putting  more  zones  in  the 
model  (here  N=250)  but  this  puts  the  eigenvalue  routine  at  a  loss  for 
solutions.  The  envelope  mass  can  be  decreased,  but  at  the  expense  of  the 
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accuracy  in  the  pulsation  modes.  The  basic  features  of  the  thermal  modes 
are  well  represented  and  accurately  calculated.  Interior  of  the 
transition  zone,  the  adiabatic  density  fluctatiuons  will  overwhelm  the 
small  luminosity  changes  due  to  the  entropy  variations. 

The  optimal  numeric  solution  to  an  eigenvalue  problem  is  found  when 
the  matrix  elements  are  approximately  equal  in  magnitude.  The 
transformation  of  the  mechanical  variable  introduced  by  Castor  (1971), 
6r-»-  XO  =  /DM2  6r  ,  does  this  for  the  adiabatic  wave  equation.  In  a  model 
where  the  shell  mass  is  a  geometric  progression  increasing  inward,  this 
transformation  is  equivalent  to  each  zone  having  the  same  acoustic 
traversal  time.  The  sound  speed  is  roughly  proportional  to  the  square 
root  of  the  temperature  and  the  thickness  of  the  shell  is  roughly 
proportional  to  the  mass;  so,  the  integrand  in  3-7  has  equal 
contributions  from  each  zone. 

This  transformation  has  the  unfortunate  side  effect  of  decreasing 
the  numeric  niceness  of  the  thermal  operator,  D.  The  thermal  diffusion 
velocity  (4-16)  decreases  as  the  sound  speed  increases.  To  optimize  the 
matrix  D,  the  mass  of  each  zone  should  decrease  inward,  the  opposite 
direction  required  of  A.  This  tradeoff  between  the  adiabatic  and  thermal 
operators  will  be  important  only  if  the  mode  is  dominated  by  thermal 
effects.  For  quasi-adiabatic  modes,  the  interior  is  dominated  by  the 
adiabatic  density  fluctuations  and  the  inaccuracy  in  the  thermal  modes 
is  unimportant.  We  are  working  on  ways  of  circumventing  this  difficulty 
so  that  both  operators  are  represented  as  accurately  as  possible. 
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Table  5-1 
BW  Vulpeculae 


M=11.0   M 

sun 

log(L/L   )=4.23 
sun 

T  ..=  25120  K 
eff 

nQ=  0^193  ,  t1q=-1.9  *  10~^ 

n^=  0^147  ,  Ti^=-2.3  *  lO"^ 

n^=   0^118  ,  Ti2=-1.5  *  10"^ 

Composition  X=0.70,  Y=0.28,  Z=0.02 


Envelope  mass  =  5.02  M 

sun 


73 


Table   5-2 
DL   Cassiopeiae 
M=5.69   M 


log(L/L  )=3.57 

sun 

T3ff=   5848   K 

nQ=   7^87    ,    tIq  =   6.8   *    lO"^ 

n^=    5^64    ,    r|^=    1.9   *    10~^ 

n^=   4^28    ,    11^=   6.2   *    lO"^ 

Composition  X=0.70,  Y=0,28,  Z=0.02 


Envelope  mass  =  2.76  M 

'^  sun 


CHAPTER  VI 
DISCUSSION 

To  understand  the  utility  of  the  thermal  eigenvectors,  we  shall 
look  at  various  forms  that  the  linearized  pulsation  equations  (3-5  and 
3-6)  acquire  when  the  ^  's  are  used  to  expand  them  in  the  ratio  of 
thermal  time  to  pulsation  period.  Eliminating  6S  one  can  rewrite  3-5  and 
3-6  as  an  integro-dif f erential  equation: 


0)2  6r  =  (  A  +  B  —^ —  C  )  6r.  (6-1) 

lU)  —  u 


we  have  an  eigenvalue  problem  for  ur  and  6r.  Eliminating  6r  ,  3-5  and  3- 
6  can  be  written  in  a  form  suitable  for  studying  the  secular  behavior  of 
a  star: 


ia)6S  =  (D  +  C — B)  6S.  (6-la) 

/,l2   _   A 


The  ioj  -  D  term  is  a  resolvent  that  can  be  expanded  in  the  complete  set 
of  eigenvectors  of  D.  Mopting  the  Dirac  (or  bra-ket,  see  Davydov  1965, 
ch.  5)  notation  of  quantum  mechanics,  this  set  is   {  ~o   »|  C   >  }  » 
and  equation  6-1  becomes 
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0)2  6r  =  (  A  +  ^^  B   .^"  ^  /  C  )  6r  .  (6-2) 


Next,  the  mechanical  eigenvector  is  expanded  in  the  set  of  eigenvalues 
and  eigenvectors  of  A,  {  w?  , |  y .  >  } .  We  would  like  to  examine  the  non- 
adiabatic  effects  on  a  mode  that  is  primarily  adiabatic,  and  this 
expansion  will  be  written  as 


6r  H  I  y  >  =  I  y.  >  +  I  a.  I  y.  > 


and 


0)  =  (0.  +  Au 


where 


(0.^   y.  >  =  A   y.  >  . 
1   '   1        '  -^1 


Inserting   this   expansion   and   multiplying   eq.    6-2   by   <   y    |    give 


(1    +   ^^jmI'^jI^    )    0)2   =   <   y.|A|    y.    >   + 


<  y.     |B|    l     X   I      |C|    y.> 
y  1    '     '      n  n    '     '    -^  1 

^n  ici)       +        a 

n 


.    <  y.     |B|    I     X   E,      |C|    y.    > 

J     Tl  J  10)      +  o 


76 


<y.  |B|  l^Xl^    |C|  y  > 

: T a.  )  + 

ICO   +   a         J 
n 


^k  ^j  ^   \ 


: ^ a.  . 

loj  +   a         J 


(6-3) 


This  expression  is  rather  cumbersome,  but  following  the  normal  quantum 
mechanical  perturbation  scheme,  we  evaluate  the  change  in  the 
eigenfrequency  using  the  initial  eigenvector.  Choosing  the  adiabatic 
mode  as  the  starting  point,  the  quasi-adiabatic  correction  can  be 
written  as 


<  y,  |B|  I     X  ^  |C|  y  > 

c,2  _  ^2  ^   I     i .  ^       ^-         (6-4) 

X    ^n         10)   +   a 


which  is  identical  to  equation  (26)  in  Castor  (1971).  The  normal  quasi- 
adiabatic  approximation  is  a  /w.  <  1  ,  for  all  possible  n,  reducing  6-4 
to 

2  _  ^2  =1    <  „_  iBCI  y.  >  (6-5) 


OJ   -  (Ji)' 


If  =  -7-^  <  y.  |BC|  y.  > 

1    10).      -^1  '    '   1 


as  y   I  C  X  C   1=1. 
^n  '   n     n  ' 


This  expression  for  the  quasi-adiabatic  correction  can  be  further 
analyzed  if  Ao)/oj  <  1.  Then  o)^  «  o)?  +  2o).Ao)  and  eq.  6-5  can  be  rewritten 


77 


Aw  =  -7-'-  <  y.  |BC|  y.  >, 
wr 


a  result  identical  to  the  y.  (equation  8,  in  Buchler  1983).  Here,  Aw  is 
purely  imaginary,  as  all  quantities  on  the  right  are  real,  and  generally 
shows  that  the  star  is  stable  (lAco  <  0)  as  the  diagonal  elements  of  B 
are  positive  and  the  diagonal  elements  of  C  are  negative. 

This  recapitulation  of  the  quasi-adiabatic  analysis  serves  to  point 
out  that  by  knowing  the  structure  of  the  thermal  modes,  a  better 
approximation  can  be  made.  Partitioning  the  thermal  modes  into  two  sets. 


and 


h   =   t  o^,    l^,    a=  1,  N^  :  0^0,.  <  1 


S,  =   I  a  ,  C  ,  n  =  N  +  1,  N  :  w,/a  <  1  }, 


eq.  6-2  can  be  rewritten  as  an  eigenvalue  problem  assuming  an  adiabatic 
starting  point  only  to  determine  N„ 

u)2  I  y  >  =  a|  y  >  + 

^  ,     io)   +   0 
n=l  n 

N     B|  I     X  &   |C|  y  > 

y     :-2— 2 .  (6-6) 

M  ^1   1(1)  +  a 
n=N^+l  n 

The  second  term  on  the  right  is  the  quasi-adiabatic  correction,  now  with 
an  appropriate  cutoff.  However,  the  cutoff  is  not  in  the  spatial 
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variable,  as  the  integrals  implied  by  the  notation  extend  over  the 
entire  mass  of  the  star.  The  correct  cutoff  is  to  stop  the  sum  at  the 
thermal  eigenvalue  whose  magnitude  is  just  below  the  pulsation 
frequency,  as  the  sum  over  the  set  Si  is  doing.  For  this  reason,  the 
second  term  should  be  more  correct  in  the  calculation  of  this  part  of 
the  stability  coefficient. 

The  third  term  produces  a  correction  in  the  lowest  order,  and 
therefore  cannot  be  treated  as  a  perturbation.  To  see  this  effect,  the 
denominators  of  the  sums  are  expanded  in  a  binomial  series  and  eq.  6-6 
becomes 

0)2  I  y   >  =  a|  y   >  + 


1  h      «l  ^n><  ^n  1^1  y> 


10)  ^ 


1  +  a  /io) 
n=l  n 

N    B|  C^  X  5^  |C|  y  > 

^  ^  ^^   a  (1  +  io)/a  )  ^^-^^ 

n=N  +1   n  n 


0)2  I  y  >  =  A|  y  >  + 


^  f     BU  XI  -^  X  C   |C|  y  >  + 
10)  ''  ,   '   n       10)     n  '  ' 


N, 

I 
n=l 

I  B|  C^Xl  -iii^X  C^  Id  y  >.     (6-8) 

n=N^+l  n 

The  first  term  is  the  quasi-adiabatic  correction  already  discussed 
and  the  second  is  the  effect  of  the  outer  layers  on  the  pulsation.  The 
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expansion  is  in  (x^l  a    ,  and  in  the  limit  w/o  — >0,  or  very  short  thermal 
'^  n  n 

times;  this  sum  becomes 

N  1 

y     Bl  5  >-^    I      |C|  y  >. 

n=N^+l         n 

This  term  represents  a  zero-order  quantity  and  should  not  be  neglected 
when  the  Initial  eigenvalue  is  calculated.  If  the  thermal  modes  had  a 
theta-function  behavior  at  the  transition  zone,  the  set  S2  would  form  a 
complete  set  in  and  of  itself  for  describing  the  thermal  response  of  the 
outer  layers  of  the  star  on  time  scales  shorter  than  o)   .  The  actual 
eigenvectors  have  a  fairly  sharp  fall  at  this  point;  so,  a  first  try  at 
including  these  terms  is  to  write  the  sum  over  states  as  D'    ,  where  D' 
is  the  outer  [(N  -  N  )*(N  -  N  )]  submatrix  of  D  that  includes  only  those 
zones  with  I  >  N   .  Substituting  D'  for  the  second  sum  over  states  gives 

0)2  I  y  >  =  a|  y  >  + 


1  \ 
-J—Y^\l     ><   I      |C|  y  >  - 
i  CO  ^    '   n     ^n  '  '  -^ 
n=l 


B|(D')"^|C|  y  >  +  ia)B|(D')   |C|  y  >.    (6-9) 

This  equation  can  be  expanded  in  the  ratio  of  the  proper  time 
scales;  the  second  terra  on  the  right  a  =  o/u„  and  the  last  two 
in  a  =   oj„/a  .  Writing  a=a  =a  to  give  the  simplest  expression  and 
arranging  the  equations  in  a  hierachy  of  terms  by  order  in  a  give 
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oo2  I  y^  >  =  (  A  -  B(D')  ^C  )|  Yq  >  (6-lOa) 


0)2  I  y^  >  +  20)^0)^  y^  >  =  (  A  -  B(D')  ^C)|  Y^  >  + 


N 

^  I     B\    I     X  C   |C|  y-  >  + 
i(i)_  ^   '   n     n  '  '  ■'0 
0  n=l 


IWqBCD')  ^C|  Yq  >;  (6-lOb) 

where 

I  y  >  =  I  Yq  >  +  a  I  y^  >  +  a2  I  y^  >  +  ... 
and 

0)  =  u)_  +  ao)  +  a2co  +  O(a^), 


To  eliminate  secular  terms  in  the  equation  for  |  Yi  >,   the  scalar 
product  <  Yq  I  y-j^  >  is  set  to  zero.  For  convenience  we  shall  set 
<yQ|yQ>=l.  This  gives  the  initial  eigenvalue  and  the  first  order 
stability  coefficient  as 


0)2  =  <  y^  |A  -  B(D')"^C|  Yq  >  (6-lla) 


-1=-;^    r<yo  |B|  l^x  l^  |c|  Yo>  + 

liiit.      n=l 

V2<  Yq  |B(D')~^C|  Yq>.  (6-llb) 

At  any  point  the  matrix  D'  can  be  replaced  by  the  sum  over  states  that 
it  approximates 


ol     =  <  Yq    \iA+  I  bU^  >  -^<  ^^  |C)|  Yq  >,  (6-12a) 

n=N„+l  n 

T 


2oj^   n=l 

n=N„+l  a^ 

T  n 

Equations  6-1  la  and  6-1 lb  agree  with  equations  25  and  39  respectively  in 

Buchler  and  Regev  (1982b),  who  derived  them  from  physical  arguments 

regarding  the  constant  luminosity  of  the  outer  layers.  It  should  be 

noted  that  the  assumption  of  ; =  0  is  not  valid  over  the  entire 

dm 

section  of  the  Cepheid  model  above  the  transition  zone.  The  sum  over 
states  is  valid,  and  comparing  the  thermal  mode  6L/L  to  that  of  the  LNA 
analysis  (Figure  5-16),  we  see  that  the  functional  dependence  of  the 
LNA  6L/L  is  reproduced  by  the  thermal  modes. 

Using  the  infinitely  sharp  cutoff  limit,  the  eigenvectors  of  6-lOa 
can  be  calculated  in  a  straight  forward  way.  Defining  the  error  in  Iyq  > 
as 


E  =  {  Ii  I  yd)  -  6r(I)  |2f/2 


where  6r  is  the  real  part  of  the  mechanical  eigenvector  given  by  the  LNA 
stability  analysis  and  both  eigenvectors  are  normalized  to  unity  in 
their  respective  scalar  products.  The  error  in  the  eigenf requency  is 
defined  as 
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w    '  0    LNA'   LNA 


These  quantities  are  plotted  in  Figure  6-1  against  the  location  of  the 

transition  zone,  loe(  T      )•  A  sharp  minimum  is  seen  at 
'   °        trans  ^ 

log(T      )  =  5.23,  corresponding  to  the  location  of  the  actual 
'^  trans         ' 

transition  zone  as  predicted  by  the  luminosity  eigenvector  plotted  in 

Figure  5-6.  In  this  model,  the  change  is  rather  abrupt,  falling  a  factor 

of  4  when  log(  T      )  changes  by  -0.13  or  +0.04  .  The  change  w_  is  not 
*   trans  '     »    •'  °        0 

as  dramatic,  reflecting  the  "minimum"  principle  that  a  quasi-adiabatic 
oscillation  will  follow. 

When  the  same  quantities  are  plotted  for  the  Cepheid  (Figure  6-2), 
it  is  obvious  that  the  agreement  is  not  much  better  than  using  the 
adiabatic  eigenvectors.  The  error  is  reduced  at  a  point  where  the 
transition  zone  is  located,  but  not  by  a  factor  of  10  as  in  BW  Vul.  This 
results  from  the  more  complicated  behavior  of  6L/L  in  the  Cepheid  (see 
Figure  5-16).  The  thermal  modes  have  a  large  peak  in  the  ionization 
zones  due  to  the  integrating  factor's  temperature  dependence  but  do  not 
have  the  rapid  decrease  inward  of  those  in  BW  Vul.  Therefore,  the 
approximation  of  replacing  the  sum  over  states  by  the  submatrix  D'  is 
not  a  good  one  for  this  class  of  star.  As  the  sum  over  states  contains 
the  necessary  information  about  the  luminosity  variations  in  the 
exterior,  including  the  integrals  as  in  6-12a  should  give  much  closer 
agreement.  This  problem  is  discussed  in  Pesnell  and  Buchler  (1983). 

Comparing  the  initial  eigensystem  6-1 2a  to  the  adiabatic  wave 
equation  produces  three  differences: 
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1)  There  is,  at  present,  no  guarantee  of  either  the  reality  or 
positivity  of  the  eigenvalues  of  either  6-lla  or  6-12a. 

2)  The  eigenvectors  <  y^  |  and  |  y„  >  are  not  necessarily 
related  by  a  transposition,  much  less  being  equal.  The  matrix 
A  -  B(D')   C  is  not  symmetric  but  is  of  an  upper  Hessenberg 
form  whose  only  advantage  is  being  diagonally  dominant.  The 
eigenvectors  come  as  right  and  left  hand  vectors  and  must  be 
calculated  separately. 

3)  There  is  no  proven  relationship  between  the  number  of  nodes 
in  I  y^  >  and  the  numeric  position  of  oo   in  an  increasing 
ordering  of  the  spectrum  of  A  -  B(D')   C  . 

However,  for  stars  which  are  predominantly  adiabatic,  as  most 
pulsators  are,  the  eigenvalues  are  close  to  the  adiabatic  frequencies, 
the  second  part  of  the  operator  having  little  or  no  effect  on  A.  The 
transposed  eigenvector  has  a  serious  numeric  problem.  In  the  region  of 
the  transition  zone,  the  transposed  matrix  has  a  structure  such  that  a 
discontinuity  appears  in  <  yQ  |.  The  use  of  eq.  6-12a  instead  of  6-lla 
will  spread  this  feature  over  several  zones,  removing  the  discontinuity. 
When  the  transposed  eigenvector  is  calculated  using  the  LNA  analysis, 
the  feature  is  present  but  is  much  smoother.  A  change  in  the  method  of 
calculation  is  indicated,  not  a  revision  of  6-1 2a. 

A  final  observation  regarding  6-1 2a  concerns  the  limit  N  — >0.  As 
a  normal  star  is  almost  adiabatic  (  N  — >N) ,  the  opposite  limit  must 
describe  an  abnormal  star.  For  the  condition  N  — >0  to  be  satisfied, 
the  non-adiabatic  regions  of  the  star  must  be  extensive,  covering  all  of 
the  star  that  is  oscillating.  These  conditions  may  be  found  in  the  R 
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Corona  Borealis  stars  or  other  highly  evolved  stars  where  the  envelopes 
are  tenuous  and  have  low  gas  to  total  pressure  ratios.  In  these  stars, 
trying  to  recover  the  adiabatic  limit  in  one  obvious  way  will  lead  to 
some  rather  strange  behavior  of  w^. 

Introducing  an  arbitrary  scale  factor,  X    ,  into  equations  3-5  and 
3-6,  which  represents  the  ratio  of  the  pulsation  period  to  some  overall 
thermal  time,  they  become 

u)26r  =  A'6r  +  B»6S  (6-13a) 

iu)6S  =  \(C»6r  +  D'6S).  (6-13b) 

For  most  stars,  when  \-K)   the  adiabatic  limit  is  recovered  as  6S  is 
forced  to  0.  However,  in  the  integro-dif ferential  form,  splitting  the 
thermal  modes  as  before , 

(i)2  I  y  >  =  A|  y  >  + 


X        A     «l  ^n><^n  \^\y>  ^ 


1   +     Xa   /ico 
n=l  n 

N    B|  IX    I      |C|  y  > 

X  I         —. 2_,j .         (6-14) 

n=N^+l  n 

The  value  of  N  may  depend  on  X   for  most  stars ,  but  for  the  radiation 
dominated  atmospheres,  the  thermal  response  time  is  very  short  for  all 
modes  whose  cutoff  in  amplitude  is  above  the  normal  transition  zone. 
Using  the  ideal  gas  as  before,  the  atmosphere  should  have  a  low  density 
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gradient  and  opacity  gradient  so  that  the  diffusion  length  (4-15)  is 
large.  If  the  outer  modes  have  a  /co  >  1  for  all  n  >  N„  ,  then  N  is 
almost  constant  during  this  discussion. 

In  the  limit  of  X   — >0  (  N  constant),  eq.  6-14  is  expanded  using 
the  binomial  theorem  as  before. 


0)2  I  y  >  =  A|  y  >  + 

N  1 

T   B|  C  >  -^  <  ^   |C|  y  >  +  0(\)  ,   (6-15) 
n=N^+l        n 

the  B(D')   C  term  does  not  disappear.  For  highly  non-adiabatic 
envelopes,  this  term  can  dominate  and  change  the  spectrum  so  that  it  no 
longer  reduces  to  the  adiabatic  case.  This  would,  in  part,  explain  the 
modes  that  have  been  named  "strange"  modes  (King  1980;  Saio  and  Wheeler 
1982). 

This  explanation  of  the  strange  modes  does  not  require  that  a  mode 
be  near  (in  some  sense)  to  an  adiabatic  mode.  There  are  no  theorems  that 
prove  that  all  of  the  eigenvalues  of  6-1 2a  are  real,  and  starting  from 
an  adiabatic  mode,  the  non-adiabatic  corrections  could  move  the 
frequency  far  away  from  the  initial  value.  It  is  interesting  that 
if  N  can  be  considered  constant,  the  adiabatic  limit  is  never  recovered 
in  some  models.  This  property  is  just  what  a  strange  mode  shows. 

Ilodes  that  possess  the  property  of  N  being  independent  of  X.  shall 
be  called  "sudden"  modes.  They  are  characterized  by  a  condition  of 
thermal  balance  existing  over  a  large  fraction  of  the  mechanical 
eigenvector's  amplitude.  The  quasi-adiabatic  stability  coefficient  would 
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be  a  poor  approximation  for  these  modes,  as  the  second  term  in  6-12a 
will  dominate  co. .  We  must  stress  that  the  condition  of  thermal  balance 
does  not  imply  that  the  luminosity  has  a  zero  derivative  in  these  modes, 
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CONCLUS IONS 

We  have  shown  that  the  thermal  structure  of  stellar  pulsations  can 
be  represented  by  the  eigenfunctions  of  the  diffusion  equation.  As  this 
equation  is  of  second  order  in  the  spatial  derivative,  it  can  be  mapped 
onto  a  Sturm-Liouville  operator,  with  the  concomitant  theorems  regarding 
the  reality  of  the  eigenvalues  and  the  form  of  the  eigenfunctions.  This 
transformed  operator  defines  a  thermal  time  scale  through  a  new 
independent  variable  with  the  units  of  (time)'2  .  This  time  scale  agrees 
with  that  of  earlier  authors  in  the  interior  and  predicts  the  same 
location  for  the  transition  zone  in  a  star. 

The  thermal  modes  that  arise  from  the  eigenvalue  problem  can  be 
used  to  analyze  the  thermal  effects  in  a  first  order  perturbation 
expansion.  The  traditional  adiabatic  wave  operator  is  replaced  by  an 
operator  that  includes  the  knowledge  that  the  outer  layers  are  in 
thermal  balance ,  and  that  this  does  not  imply  that  ^  ^  =  0  in  these 
layers. 

Partial  spectra  of  the  thermal  modes  are  presented  for  two 
realistic  stellar  models,  along  with  diagrams  of  the  eigenvectors. 
Mathematically  consistent  expansions  of  the  pulsation  equations  are 
presented  that  give  the  quasi-adiabatic  stability  coefficient  and  the 
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analogous  expression  when  the  outer  layers  are  correctly  accounted  for. 
It  is  shown  that  a  new  type  of  mode,  the  "sudden",  mode  can  be  defined 
using  these  expressions. 

The  driving  mechanism  in  Cepheids  has  features  in  the  thermal  modes 
that  do  not  show  up  in  a  star  of  the  Beta  Cephei  class.  This  casts  doubt 
on  the  viability  of  ionization  driving  being  responsible  for  these 
stars'  variability. 

In  an  aesthetic  vein,  a  way  of  displaying  the  entropy  or 
temperature  variations  for  stars  that  normally  have  a  "spike"  in  the 
hydrogen  ionization  zone  is  given,  showing  much  of  the  detail  in  this 
region  without  the  introduction  of  an  arbitrary  cutoff.  This 
transformation  also  increases  the  numeric  accuracy  of  solutions  in  the 
linear  non-adiabatic  pulsation  stability  analysis  by  limiting  the 
variation  of  the  elements  of  the  coupling  matrices. 


APPENDIX  A 
THE  ANALYTIC  EQUATION  OF  STATE 

For  the  purposes  of  these  calculations,  an  equation  of  state  has 
been  developed,  based  in  part  on  a  routine  supplied  by  Arthur  Cox  (Cox, 
A.  N.  1980b).  In  this  appendix  the  various  formulae  and  constants  used 
are  discussed. 

As  our  analysis  requires  the  entropy  of  an  ionizing  gas,  we  first 
calculate  the  free  energy  ( F)  as  a  function  of  the  specific  volume  (v), 
temperature  (T) ,  and  composition  (X,Y,Z).  All  thermodynamic  functions 
can  be  evaluated  by  derivatives  of  F  (Zel'dovich  and  Razier  1966,  ch. 
Ill): 

entropy:  S  =  -  (y^ 

'  v,N 

pressure:  P  =  -\^j^^ 

0/5  (F/T)\ 
internal  energy:  E  =  -  T^l — :—z — |   . 

\      A.N 

In  order  to  calculate  F,  the  approximation  of  an  ideal,  ionizing 
gas  is  adopted  with  the  characteristics  of  the  five  species  listed  in 
Table  A-1.  Values  for  the  mass  weighting  of  the  elements  has  been 
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adapted  from  Cox  and  Tabor  (1976),  and  the  atomic  masses  and  ionization 
potentials  are  from  Novotny  (1973).  For  completeness,  a  term  is  added  to 
include  the  effects  of  black  body  radiation.  This  is  assumed  to  be  that 
appropriate  to  isotropic  radiation  so  that 


rad     3 


The  ionized  fraction  of  each  element  is  determined  by  minimizing 
the  free  energy  at  the  given  temperature  and  specific  volume.  This 

procedure  leads  to  a  Saha  equation  that  is  solved  using  a  Newton  I-tethod 

12 
to  a  relative  accuracy  of  1  part  in  10   for  the  electron  density.  When 

the  electron  density  is  known,  the  reciprocal  of  the  mean  molecular 

weight  (in  moles/gram)  is 

n  =  1/ti  =  X/1. 00797  +  Y/4.0026  +  Z/21. 02447  +  n^. 

The  relative  abundance  of  each  ionization  stage,  n^ ,  is  calculated  from 
the  ratios  of  the  appropriate  partition  functions,  and  the  free  energy 
is  written  as  a  sum  over  these  species. 

F  =  -  RT  {  y.  n,ln(m^^^u./n.)  +  n[  -2.5048  +  ln(vT-^^^)  ]  } 
'■  ''1  i    1   1  1     '- 


-  -  vT'* 
3 


where  the  constant  term  is ; 


95 


1  +1-  ln(2ran  k)  -  ln(N.)  -3  ln(h), 
2        p  A 


m.  is  the  species  gram  molecular  weight,  and  u^  is  the  temperature 
independent  partition  function  that  is  assumed  to  the  ground  state 

degeneracy. 

2 
The  expressions  for  the  entropy  (ergs/gram/K) ,  P  (dynes/cm  ),  and  E 

(ergs/gram)  are 

S  =  r{1^   n.ln(m^^^u./n^)  +  n[  -1.0048  +  ln(vT^^^)  ]  }  +  -^  vT^ 


a  4 
P  =  nRT/v  +  ^  T 


E  =  Y  nRT  +  avT**  +  l^   n^  x^ 


Derivatives  of  the  pressure  and  the  internal  energy  are  found  by 
analytic  differentiation  of  the  above  formulae. 

The  subroutine  NEDS  performs  the  calculations  required  to  find  the 
thermodynamic  quantities  as  outlined  above.  A  listing  of  this  routine 
follows. 
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TABLE    A-1 


PHYSICAL   PROPERTIES    OF  THE    COMPONENTS    OF  THE   GAS 


f/ 

species 

atomic  weight 

u 

X 

(gr/raole) 

(eV) 

1 

Neutral  Hydrogen 

1.00797 

2 

13.595 

2 

Ionized  Hydrogen 

1.00797 

1 



3 

Neutral  Helium 

4.0026 

2 

24.581 

4 

Ionized  Helium 
(singly) 

4.0026 

4 

54.403 

5 

Ionized  Helium 
(doubly) 

4.0026 

2 



6 

Ionized  metal 

24.969 

1 

5.524 

7 

Ionizing  metal 
(single  stage) 

45.807 

1 

7.900 

8 

Neutral  metal 

17.969 

- 



9 

Electrons 

5.486(-4) 

2 
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APPENDIX  B 
THE  INITIAL  MODEL  INTEGRATOR 


For  the  present  calculations,  an  initial  envelope  model  gives  an 
adequate  description  of  the  variation  of  the  physical  parameters  as  a 
function  of  position  in  a  star.  To  construct  an  initial  model  we  must 
solve  the  equations  governing  stellar  structure  with  time  derivatives 
set  equal  to  zero.  With  the  appropriate  boundary  conditions,  a  solution 
can  be  found  by  integrating  inward  from  the  surface  until  the  desired 
envelope  mass  is  reached.  In  this  appendix,  the  equations,  boundary 
conditions  and  integration  techniques  are  discussed. 

A  spherically  symmetric  model  in  hydrostatic  and  thermal  balance  is 
governed  by  the  equations 


d  r^ 

conservation  of  mass    — —  =  — ; — rr  (B-1) 

d  m  4ti/3 


d   P  G  M  ,^   „, 

conservation  of   momentum       — —  = r  K'i~^) 

dm      ,4 
4ur 


conservation  of  energy  -: —  =  0  (B-3) 

dm 


4 

T/  X     //   2.2  ac  d  T 

energy  transport  LCm)  =  -  (4-nx  ;  — :; ;; — 

^■'  '^  3k  d  m 
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(B-4) 


104 
where  they  have  been  written  in  a  Lagrangian  form  that  uses  the  interior 
mass  as  the  independent  variable.  The  energy  transport  is  assumed  to  be 
purely  radiative  and  the  opacity  is  a  fit  due  to  Stellingwerf  ( 1975a, b). 

We  divide  the  model  into  N  mass  shells,  or  zones,  with  each 
variable  assigned  to  the  zone  center  or  zone  interface  as  shown  in 
Figure  B-1.  The  mass  of  shell  I  (zone-centered  mass)  is  called  DMl(I), 
and  the  mass  associated  with  interface  I  is 

DM2(I+1)  =  V2  (DM1(I)+DM1(I+1)). 

A  mass  derivative  of  a  zone-centered  quantity  (pressure  etc.)  is 
interface  centered  and  vice-versa.  All  quantities  are  seen  to  be 
correctly  centered  with  the  exception  of  the  opacity  in  equation  B-8. 
The  presence  of  the  zone-centered  opacity  in  the  equation  for  the 

interface-centered  luminosity  necessitates  the  use  of  special 

4 

interpolation  schemes  to  correctly  center  the  —  -^ —   term.  Several 

*^  ^  K  d  m 

authors  have  published  formulae  for  this  purpose  and  we  have  chosen  that 
of  Stellingwerf  (1975a)  for  our  calculations.  The  specific  form  of  the 
difference  equations  is  such  that  they  are  compatible  with  the  nonlinear 
hydrodynamic  code  developed  at  Los  Alamos  National  Laboratory  (King  et 
al.  1964;  Cox,  Brownlee ,  and  Eilers  1966).  This  work  forms  the  basis  for 
further  work  in  nonlinear  behavior  of  models  and  this  compatibility 
enables  us  to  compare  our  results  with  those  of  fully  nonlinear  models. 
With  this  convention  for  the  masses  and  derivatives,  equations  B-1 
through  B-5  can  be  written  as 
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Rd+l)-^  -  R(I)^  _  v(I)  .„_.. 

DM1  (I)        4ti/3  •  ^      ^ 


P(I+1)  -  P(I)  ^    G  M(I+1) 
DM2(I+1)    -  "4,  R(i+i)4 


(B-6) 


and 


L(I-M)  -  L(I)  ^  _ 

DM1  (I) 

4 

L(I+1)  =  -(4tiR(I+1)^)^  -^  <-^>.  (B-8) 

J     K  a  m 


The  boundary  conditions  at  the  surface  are  P=Pq  and  T=TQ=0.841Tg^£ , 
where  Pq  is  the  pressure  due  to  an  external  atmosphere  and  T  ^^  is  the 
effective  temperature  of  the  model.  As  these  conditions  are  valid  at  a 
point  beyond  the  actual  model  surface,  they  must  be  incorporated  into 
the  differenced  equations  by  giving  some  elements  of  B-6  and  B-8  special 
forms. 

In  the  momentum  equation  (B-6),  the  pressure  derivative  is 
rewritten  as 


P(W-l)  -  P(N)  _      P(N)  .   Qs 

DM2(Nfl)   -  y  • 

When  f .=  1   the  outer  pressure  (P(I«-1)=Pq)  is  zero,  but  there  is  a 

better  approximation  for  f   (see  King  et  al.  1964).  If  one  assumes  that 

an  infinite  number  of  zones  extends  outward  from  the  surface  with  each 

zone  mass  a  constant  fraction  of  the  one  interior,  f  .   can  be  written  as 

A 


f  =    "  ^  \  (B-10) 

A    a  +  1 
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where  a  =  DM1(N)/DM1(N-1)  is  the  constant  mass  ratio.  Values  for  f  are 
normally  0.1-0.2  .  This  modification  moves  the  zero  pressure  boundary 
condition  to  the  outermost  zone  as  desired.  For  a  consistent  treatment, 
radiation  pressure  is  added  to  this  "atmospheric"  contribution. 

When  the  boundary  condition  T  =  Tq  is  inserted  into  the  luminosity 
equation  (B-7),  an  analogous  expression  can  be  derived  for  the 
temperature  of  the  outermost  mass  layer.  Starting  with  equation  B-4 ,  and 
assuming  the  luminosity  and  radius  are  constant,  a  grey  atmosphere 
relationship  for  T(N)  is  found: 


t'*(N)  =  T^^^  I  (  T  +  2/3  )  (B-11) 


where 


X  =  Jk  dm/4-rtR(Nfl)2  (B-12) 


Given  T    ,  DMl(N),  R(N+1),  and  the  opacity,  T(N)  can  be  evaluated  so 
that  the  boundary  condition  is  satisfied. 

Along  with  the  assumed  form  for  the  outermost  mass  shell,  there  are 
two  additional  relationships  that  need  to  satisfied  to  give  a  physical 
model.  One  of  these  is  that  the  radius  evaluated  at  the  effective 
temperature  be  the  photospheric  radius  given  by  the  black  body  formula 


L  =  4ti0  R  J  ^   /..  .  (B-13) 

photo  eff 


This  is  satisfied  by  introducing  an  auxiliary  expression 
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R(W-l)  =   R  .    (  1  +  P  )  (B-14) 

photo 


and  varying  P  so  that 


R(  T(m)  =  T^^p   =  Rp^^^^ 


The  second  condition  requires  the  optical  depth  of  the  outermost 

-3         -2 
shell  be  at  a  chosen  value,  x„  ,  where  10  <  x^  <  10   .  This  is 

achieved  by  varying  the  outer  zone  mass  until  x(  eq.  B-12  )  =  x^  .  Both 

of  conditions  are  satisfied  to  a  relative  accuracy  of  one  part  in  one 

million. 

Given  eq.  B-5  through  B-8  and  boundary  conditions  B-9 ,  B-11,  Br-12, 

and  B-14,  we  have  a  complete  set  of  equations  to  describe  the  envelope. 

To  integrate  this  set,  the  mass,  luminosity,  effective  temperature,  and 

composition  are  read  in,  guesses  for  DMl(N),  f^  ,  and  P  are  made,  then 

values  for  P(N)  and  T(N)  are  found  from  B-9,  B-6  and  B-U.  This  zone  is 

automatically  in  thermal  balance  and  v(N)  is  changed  to  satisfy 

hydrostatic  balance  in  a  one  dimension  Newton  method.  The  optical  depth 

of  the  zone,  x(N)  ,  is  compared  to  x^  ,  and  if  they  do  not  agree,  the 

mass  of  the  zone  is  modified  following  another  Newton  method  until 

x(N)  =  x„  to  the  desired  accuracy. 

After  the  outer  zone  is  converged,  eq.  B-5  is  used  to  find  the 

inner  radius  R(N).  The  mass  for  zone  N-1  is  a  constant  multiple  of 

DMl(N),  and  this  constant  can  be  varied  so  that  T(N-1)/T(N)  <  y,   where 

Y  is  normally  1.06.  In  a  gas  whose  dominant  component  is  hydrogen,  the 

electron  density  changes  rapidly  during  the  ionization  of  hydrogen.  This 
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Increase  in  electrons  corresponds  to  a  rapid  increase  in  the  opacity  as 
well.  A  star  whose  surface  temperature  is  greater  than  15,000  K  will  not 
have  this  difficulty  as  the  ionization  of  helium  does  not  produce  as 
rapid  an  increase  in  the  electron  concentration.  By  limiting  the 
temperature  gradient  as  above,  the  opacity  and  opacity  derivatives  are 
not  permitted  to  vary  so  much  that  they  introduce  further  inaccuracies 
into  the  calculations. 

Once  DM1 (N-1)  is  known,  the  hydrostatic  pressure  necessary  to  keep 
DM1 (N)  in  suspension  is  given  by  eq.  B-6.  A  first  guess  for  T(N-l)  and 
v(N-l)  is  made  and  they  are  corrected  by  a  two  dimension  Newton  method. 
Convergence  of  T  and  v  to  the  desired  accuracy  ends  the  iteration  in 
this  zone.  This  process  of  guessing  (if  necessary)  the  mass  ratio, 
evaluating  the  hydrostatic  pressure  and  zeroing  the  functions 


|P   -  P.  ,   I  and  IL   t,  o~  L  ^ 
'  eos   hydro'      '  eq.B-o   star 


is  repeated  for  each  succesive  interior  zone  until  the  desired  envelope 

mass  has  been  integrated. 

A  subsidary  loop  controls  the  radius  relationship.  Using  the 

initial  guess  for  p  ,  eq.  B-5  to  B-8  are  integrated  to  a  point  where  the 

temperature  exceeds  the  effective  temperature.  By  linear  interpolation 

in  the  mass  shells,  the  radius  at   T  ^^  is  found  and  compared 
'  eft 


photo  eff     photo'     phot 

true,  p  is  changed  using  a  secant  method,  R(Nfl)  is  recalculated,  and 


to  R  .  ,   .  While  the  condition  |  R(T  ..)  -  R  H^^^l  /  ^r^hoi-n>  ^°   ^® 
ohoto  eff     photo     photo 


the  integration  restarted. 
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At  the  end  of  the  integration,  all  equations  are  satisfied  to  the 
desired  accuracy  (here  one  part  in  one  million.)  This  integration  has 
been  extremely  stable  for  all  physically  viable  models  and  no  traps  or 
inconsistencies  are  known  to  exist.  The  routine  ATMHAT  and  support 
routines  for  the  model  Duilder  are  listed  below. 
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